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SYNOPSIS 


I'he present study describes the liiieai' spatial 
stability characteristics of the developing flow in a two 
dimensional parallel plate channel and in a rigid circular 
pipe while accounting for all non-parallel effects. It also 
compares the par'allel.-flow stability characteristics of the 
developing velocity profiles obtained from the linearization 
and finite difference methods. Both symmetric and antisymme- 
tric disturbances are considered for the parallel flow stability 
characteristics of the developing flow in a chaimel while only- 
symmetric disturbances are considered when de-termining the 
non-parallel effects o’n the stability characteristics since 
they are more unstable , Por -'Ghe developing flov/ in a pipe , 
only axisymmetric disturbances are considered. 

The method of multiple scales has been used for 
studjing the stability characteristics of the non-parallel 
developing flow in these ducts. This method leads to a 
sequence of ordinary differential equations, the lov/est (zeroth) 
order of which is the familiar Orr-Sommerfeld aquation 
describing the parallel flo-w stability problem. The present 
study restricts itself to the first order non-parallel 
effects only since the higher order effects are expected to 
be negligible owing to tho large Reynolds number involved 
for both the flows. Instead of integrating Oi-r-Sommerfeld 
equation for the parallel flow stability characteristics, 



this analysis integrates the continuity and momentum equations 
separately for the disturbance Telocity components and pressure 
since this method is found to be both economical and more 
accurate. The fourth order Sunge-Kutta method is used along 
with a selective application of the G-ram-Schmidt orthonormal - 
ization technique in order to remove the parasitic growth of 
error during integration. Muller's iteration has been used 
for finding the eigenvalues. 

Results for parallel flov/ stability are summarized 

below ; 

(i) The developing flow in a channel is more unstable and 
has lower critical Reynolds numbers for symmetric 
disturbances as compared to antisymmetric disturbances. 
There are, however, some eigenstates that lie close to 
the lower branch of the neutral curve and correspond to 
Reynolds numbers much greater than the critical value 
for V'Thich the antisymmetric disturbances are more 
unstable. Ror the antisymmetric disturbance, the 
critical Reynolds number R^ first decreases to about 
12100 at X cr 0.005 and then increases rapidly so that 
the fully developed flow seems to be stable to such 
disturbances. Here X = x/aRj a being half width of 
the channel (or radius of the pipe), x being the 
streamwise distance from the entry section and R being 
the Reynolds number based on a and the average velocity 
of flow. 


(ii) For the channel flow the critical Reynolds numher, ijo<piveY\umhe>r 
OLvvi frequency asad wenrenuB fe sr for symmetric disturbances 

decrease with increasing X and approa,ch asymptotically 
the corresponding Ta.lues of 3848.1, 1.0198 a-nd 0.40369 
respectively for the fully developed flow. For pipe 
flow, however, the critical frequency and wavenumber 
decrease continuously with X but the critical Reynolds 
number decreases to a value of about 11700 at X 
0.0035 and then increases monotonically to infinity. 

(iii) The stability characteristics are very sensitive to 

the velocity field description of both the flows 5 the 
velocity field obtained from the finite -difference 
method being more unstable than that obtained by the 
linearization method . In the near -entry region of 
both the flows, the critical Reynolds number for the 
velocity obtained bjr the finite-difference method is 
about half of that for the velocity field determined 
by the linearization method. 

Results for the actual non-parallel flov; stability 
are summarized below. 

(i) The actual developing flow in both geometries becomes 
unstable at lower Reynolds numbers and over a wider 
range of frequencies as compared to its parallel flow 
approximate . 

(ii) The growth rate and vmvenumber of the disturbance 

become functions of X, Y (the transverse coordinate), 



and the disturbance property for the non-parallel flow. 
Vifhile the non-parallel effects on the growth rate of 
any disturbance property are substantial j they are 
negligible for the wavenumber obtained from the parallel 
flow theory. The groTrth rate of the disturbance stream 
function is found to be maximum at Y = 0 (the channel 
centreline or the pipe axis) for all values of X, 
Reynolds number and frequency. 

(iii) Tor both the flows neutral curves at various X are 

provided for the stream function, the streamwise compo- 
nent of velocity and the energy density of the distur- 
bance. Feutral curves obtained from (X, 0), the 
growth rate of the disturbance stream function at Y = 0, 
give the minimum R and encompass the maximum unstable 
region. 

(iv) In comparison to the results for the parallel flow 

theory, the critical frequency is larger whereas E is 
smaller for all disturbance poroperties at all X for 
both the flows; however, the behaviour remains similar. 
In comparison to the value obtained for gvy (X, 0) by 
the non-parallel flo?/ theory, the parallel flow theory 
overpredicts R by 23^ for iiie channel flov/ and by 
27 % for the pipe flow at X = 0.001. Based on gip(X, 0) , 
the first instability of the developing flow in a pipe 
appears at x/a cx 33 while the parallel flow theory 
yields a value of x/a cx 37 . 



xziii 


(v) For the' pipe flow, the vs. X curve obtained on the 
basis of (Xs 0) is closest to the experimental 
data of Sarpkaya. Such a comparison for the channel 
flow is not possible since no experimental work is 
available for the same. 
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INTRODUCTION 

Following Reynolds experimental observation [l]/ 
several attempts have been made to investigate the stability 
of laminar flows to small disturbances. This stability 
analysis may be either temporal ur spatial, in the temporal 
stability analysis / disturbances which are periodic in the 
downstream direction are assumed to be applied at an initial 
instant everywhere in the fluid and are observed as time 
elapses. In the spatial stability analysis, disturbances 
which are periodic in time are imposed at a specific location 
in the fluid and are observed during their propagation down- 
stream. The flow is considered to be stable, neutrally stable, 
or unstable depending on whether these disturbances decay, 
remain unchanged, or amplify with respect to time for the 
temporal case or with respect to the downstream distance 
for the spatial case. 

1.1. Literature Survey 

The following four sub-sections critically discuss 
the available literature for the flow and its stability in 
a two dimensional parallel plate channel and in a circular 


pipe . 
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1.1.1. parallel Plate channel 

The plane Poiseullle flow has been found to exhibit 
instability at Reynolds number, based on half width of the 
channel and average velocity of flow, of about 3850 as 
obtained theoretically by Lin [2], Thomas [3], Hains [4], 
Orszag [s ], Chock and Schechter T^], and Davey [7] from the 
linear stability theory. Experiments conducted by Nishioka 
et al. 1^8] and Karnitz et al. [9] support the predicted 
results of the linear theory. However, Sherlin [10 ], 

Narayanan and Narayana [11], Patel and Head [12], and Breslin 
fl3] found the critical Reynolds number to be 3000 or less in 
their experiment probably because the disturbance level in 
the upstream flow was reasonably high. 

The temporal stability characteristics of the hydro- 
dynamically developing flow in the entrance region of the 
channel have been investigated by Hahneman et al, '][14:] and 
by Chen and Sparrow [l5] . Hahneman et al, computed the 
critical Reynolds numbers by applying the approximate formulae 
deduced by Lin [16] from the zeroth=order asymptotic solution 
of the disturbance equations. These formulae do not provide 
accurate results in the developing regicn because conditions 
postulated in their derivation are not fulfilled. Moreover, 
they used Schlichting' s [17] velocity profile for the 
mainflow which does not give accurate values of the velocity 
gradients. As is well known, such gradients play a signi- 
ficant role in the behaviour of the stability characteristics. 
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This probably explains the large scatter among the calculated 
points in the plot of critical Reynolds number versus axial 
location shown by Hahneman et al. 

Chen and Sparrow [15] used each of the three viscous 
solutions (singular, regular, and composite), obtained by 
asymptotic method, along with the inviscid solution to 
provide the critical Reynolds numbers and neutral stability 
curves by taking the basic flow velocity as (i) channel 
velocity profile, and (ii) a boundary layer velocity profile. 
As is demonstrated later, none of these solutions give 
accurate neutral stability curves. They also used two - 
numerical methods to obtain critical Reynolds numbers. While 
Nachtsheim's method [is] was found by them to be quite 
involved and it also failed at high Reynolds numbers, the 
finite difference method [3] was used for conputing only the 
critical Reynolds numbers. Therefore, the neutral stability 
curves reported by them are not accurate. Further they used 
the linearization method of Sparrow et al. [19] to calculate 
the channel velocity profile, hereafter, referred to 
as the Sparrow profile. However, it has now generally been 
accepted (cf. Section 1.1.4) that the Bodoia and Osterle 
method [20] gives better velocity field description. Never- 
theless no attempt has heretofore been made to study the 
stability of the velocity profile (hereafter referred to as 
the B-0 profile) obtained by the Bodoia and Osterle method. 

It is also well known [21-23] that the stability of parallel 
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flow is very sensitive to changes in the basic velocity 
profile. Even a little change in the velocity profile may 
change the neutral stability characteristics to a large extent 
Therefore, it becomes all the more important to study the 
stability characteristics of the developing flow in a channel 
for the B-o profile. 

Moreover, Chen and Sparrow's parallel flow approxi- 
mation of the developing flow in the channel is not correct 
since the flow is actually ncnparallel. Shen [24] has shown 
that the effect of nonparallelism of the flow is to widen the 
unstable region and to reduce the critical Reynolds number. 

In the case of a channel, one expects the effect of non- 
parallelism of the flow in the near entry region to be of the 
order of that for the boundary layer flow over a flat plate 
since the flow in this region of the channel can be treated 
as one of boundary layer type. This is expected since 
thickness of the boundary layer is much less than that of the 
central core in this region. Experiments of Schubauer and 
Skramstad [25] and Ross et al. [26] for the boundary layer 
flow over a flat plate show that the parallel flow theory 
predicts about 30% higher critical Reynolds number. The 
effects of nonparallelism of the flow should, therefore, 
necessarily be considered while studying the stability 
characteristics of the developing flow in the channel. 
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1.1.2. Rigid Circular Pipe 

The temporal stability o£ Hagen-Poiseuille flow to 
infinitesimal, axisymmetric disturbances was, perhaps, first 
analysed by sexl [27] and later by Pretsch [28], Pekeris [29] 
Corcos and Sellars [so] , Schensted [3l], and Davey and Drazin 
[ 32 ] while for the non -axisymmetric disturbances the analysis 
has been dene by Lessen, saddler and Liu [33], Burridge [34], 
Graebel [35], and salwen and Grosch [36]. For the linear 
spatial stability characteristics, theoretical analyses by 
Gill [37,38] and Davey and Drazin [32] and experimental study 
by Leite [39] for the axisymmetric disturbances, and theore- 
tical analyses by Garg and Rouleau [40] and Gill [38] and 
experimental study by Lessen et al, [ 41 ] for the non- 
axisymmetric disturbances are available. Except for the work 
of Graebel [35] and Lessen et al. [4l] , all these studies 
reveal that the Hagen-Poiseuille flow is stable to all inf ini 
tesimal disturbances. While Graebel' s analysis has errors 
[36] , Lessen et al. had too large a disturbance level for 
the linear stability results to hold. 

The temporal stability behaviour of the hydrody- 
namically developing flow in the entrance region of the pipe 
to the infinitesimal, axisymmetric disturbances was first 
studied by Tatsumi [42,43]. He used an asymptotic series 
solution and showed that instability of the flow exists in 
the entrance region of a circular tube. Tatsumi found that 
the critical Reynolds number decreases from infinity at the 
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tube inlet to a minimum at some location within the entrance 
length and then increases monotonically to infinity further 
downstream. However, Tatsumi obtained the mainflow velocity 
profile by considering the developing flow as an almost 
similar flow. This velocity field description is inferior to 
the one obtained by more recent methods (cf. Section 1.1.4). 

It has also been pointed out by Chen [44] that there is an 
error in Tatsumi' s calculation of the velocity field. Besides, 
Tatsumi used an essentially boundary layer model for the 
stability calculations in the entrance region of the pipe. 

Chen and Sparrow [is] have, shown that the boundary layer model 
does not provide accurate stability results. Therefore, 

Tatsumi 's results are not reliable. 

Huang and Chen [45,46] also analysed this problem 
on temporal basis while using Sparrow's profile, the velocity 
profile obtained from the linearization method of Sparrow et 
al. [ 19 ]. They considered both axisymmetric and non-= 
axisymmetric disturbances and found that in the near entry 
region axisymmetric disturbances are more unstable than the 
non -axisymmetric disturbances. However, their choice of 
Sparrow's profile for the mainflow is open to question as 
literature survey shows that the Hombeck method [47] is 
more accurate. No attempt has heretofore been made to study 
the stability of the velocity profile (hereafter referred 
to as the Hornbeck profile) obtained by the Hornbeck method. 
This is important in view of the fact that the critical 
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Reynolds numbers found by Huang and Chen are much larger than 
the experimental values obtained by sarpkaya [48] in the 
developing region , 

Except for an extension of Huang and Chen's work 
[ 45 ] by shen et al. [4 9], no other work is available for the 
stability characteristics of the actual (nonparallel) developing 
flow in a pipe. However, they considered only the effect of 
the radial component of velocity on the otherwise parallel 
flow temporal stability analysis and found that the minimum 
critical Reynolds number dropped to 19670 from the value of 
19900 for the parallel flow theory. They did not consider 
such nonparallel effects as the effect of streamwise variation 
of the wavenumber, eigenfunctions and growth rate, 

1.1,3. Nonparallel Effects 

The idealization of a given flow by a parallel or 
quasi-parallel flow for the purpose of studying its linear 
stability characteristics is very common. In reality, only 
fully developed, two dimensional flows are exactly parallel. 

The shear layers such as jets and wakes, all bound ary -layer 
type flows including developing flow through ducts have 
velocity components that are functions of both the streamwise 
and transverse coordinates. The parallel flow hydrodynamic 
stability theory has had only qualitative success in such 
cases. Experiments of sato and sakao [50], Ross et al. [26], 
Scotti and Corcos [Sl] , and Mattingly and Criminale [52] 
on the stability characteristics of the jet flow, boundary 
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layer flow over a flat plate, stratified shear layer and two 
dimensional wake, respectively, show systematic differences 
from the quasi-parallel flow theory. 

The difficulty in solving for the stability charac- 
teristics of a nonparallel flow lies in the evaluation of the 
eigenvalue of a set of partial differential equations. The 
basic problem is that the equations for linear stability of 
the two dimensional flows are nonseparable . To overcome this 
difficulty it is generally assumed that the nonparallel effects 
are of high order. Benney and Rosenblat [53] were perhaps 
the first to suggest that the method of slowly varying 
approximation be applied to flow stability problems, Lanchon 
and Eckhaus '[54] considered the stability of ncnparallel flows 
by using straightforward expansion about the parallel flow 
solution. However, they did not present any quantitative 
results. Barry and Ross [55], Boehman [56], Haaland [57], 
wazzan et al. [58] , Bajaj and Garg [59] and shen et al, [49] 
accounted for some nonparallel effects by retaining only the 
transverse (or radial) component of velocity. In this model, 
the governing equations are separable and they reduce to a 
modified Orr-Sommerf eld equation. 

Joseph [60] used the modified Orr-Sommerf eld 
equation as the zeroth -order term in his perturbatioi scheme 
to study the nonlinear stability of nearly parallel flows. 
Voldin [61] presented an' expansion procedure which accounts 
for all nonparallel effects. However, he assumed the lowest 
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order problem to be separable. Ling and Reynolds [6 2] also 
presented a perturbation scheme that accounts for all non=- 
parallel effects. However, they considered the temporal 
rather than the more realistic spatial stability problem, and 
expanded the stream functions of the basic flow and of the 
disturbance along with the eigenvalues in a power series about 
a given streamwise location. Their analysis is, therefore, 
valid only for a small neighbourhood of that streamwise 
location. It also neglects the vertical structure of Orr- 
Scmmerfeld solutions [63]. 

Bouthier [64,65] developed the method of multiple 
scales for steady, spatially dependent shear flows and applied 
it to the stability of the boundary layer on a flat plate. 
Other studies that consider all nonparallel effects are by 
Nayfeh, saric and Mook [66] , caster [63], Saric and Nayfeh 
[67], and Smith [68] for the Blasius boundary layer, by 
Drazin [69] for the flow in a channel whose width is a slowly 
varying function of time and streamwise coordinate, by Eagles 
and Weissman [70] for divergent channel flow, by Crighton 
and Caster [71] for developing axisymmetric jet, by Eagles 
[ 72 ] for slowly varying flow between concentric cylinders, 
and by carg and Round [73], and carg [74] for the Bickley 
jet. Sane criticism of the work by Nayfeh et al, and by 
Drazin is also available in Eagles and Weissman, For the 
present study of ncnparallel effects, we essentially use 
the approach of Eagles and Weissman. 
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1.1.4. velocity Field 

The developing flows in the channel and pipe, 
because of their practical importance, have been widely 
studied in the past and are still attracting many researchers. 
The methods used so far for calculating the velocity and 
pressure fields can be classified broadly as (i) finite 
difference solution of the boundary layer equations [20,47, 
75,76}, (ii) linearization of inertia terms [l9,77'=80], (iii) 
integral methods [81“84], and (iv) series expansions [85=87],. 
Series expansions do not provide velocity gradients with good 
accuracy [14] and one does not expect very accurate results 
from integral methods in comparison to those obtained from 
direct integration of the differential equations. The closed 
form expression for velocity and its derivatives obtained by 
Sparrow et al, [19] using linearization technique is definitely 
easier to use but the literature survey reveals that the 
velocity profiles obtained by Bodoia and Osterle [20] for the 
channel and by Hombeck T47] for the pipe are more accurate. 
Bodoia and Osterle devised a finite difference method for 
obtaining the velocity field for the developing flow in a 
channel and Hombeck extended it to the case of developing 
flow in a pipe. 

Sparrow et al. also agree that the solution 
obtained by the finite difference method is likely to be more 
accurate. The work of Crane and Burley [88], Schmidt and 
Zeldin [89] and Shah [90] clearly reveals that the Hombeck 
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method gives more accurate velocity description^ length of 
entrance region^ and pressure drop in the developing region 
of the pipe. For the channel flow van Dyke [91] considers the 
3-=0 profile to he more accurate. This view has also been 
supported by Morihara and Cheng [92] ^ Sh^h [90 ], on the basis 
of data reported in [93], concludes that the B-O profile is 
the most accurate one. Schmidt and Zeldin's [89] comparison 
of the pressure drop coefficient obtained by various methods 
supports further the view that the Bodoia and Osterle method 
is superior. 

1.2. Present Investigation 

From the above survey it is found that the stability 
analysis of the developing flow in a channel and pipe has not 
been made for the velocity profile obtained by the finite 
difference method which describes the velocity field in these 
two domains more accurately in comparison to that obtained 
by other methods. Since the stability characteristics of the 
flow are sensitive to the velocity field description [21-23], 
it is important to carry out the above analysis. Also, except 
for the work of Shen et al. [ 49 ] in which they considered only 
the effects of radial component of the mainflow velocity on 
the stability characteristics of the developing flow in a 
pipe, no study has heretofore been made for determining the 
extent to which the nonparallelism of the developing flow 
affects the stability behaviour of the flows. The above 
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literature survey reveals that these effects may be as much 
as 30%. Moreover, only the temporal stability of the deve- 
loping flow in a channel and pipe has heretofore been investi- 
gated whereas the spatial stability analysis is more realistic. 

The present work, therefore, aims to study the 
linear, spatial stability of the developing flow in a two 
dimensional channel and in a circular pipe for the velocity 
profile obtained by the finite difference method. While all 
non parallel effects on the stability characteristics of both 
the flows are studied, the infinitesimal disturbance is 
assumed to be axisymmetric for the pipe flow stability and 
symmetric for the channel flow stability. A comparative 
study of the parallel flow stability of the velocity profile 
obtained for both the flows by the finite difference and 
linearization methods has also been carried out for these 
types of disturbances. Moreover, in order to check Grohne ' s 
conclusions that only the symmetric mode is likely to lead to 
instability of the plane Poiseuille flow [94] , the velocity 
profile obtained for the channel flow by the finite difference 
method has also been investigated for its parallel flow 
stability to antisymmetric disturbances. 

For the pipe flow stability, main consideration 
is given to the wall mode since for developing flow this is 
found to be more unstable than the central mode [46] . This 
is true because the developing flow is of the boundary 
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layer type and the instability of the flow^ if it exists, 
should originate near the pipe wall as in the boundary layer 
flow. For the sake of completeness, however, the central 
modes are also studied for the parallel flow stability. 
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Chapter 2 

THEORETICAL ANALYSIS 


In this chapter stability equations which include 
all nonparallel effects have been obtained for the developing 
flow in a circular pipe and in a two dimensional channel. This 
is followed by a discussion on the appropriate boundary 
conditions. 

2.1. Development of the stability Equaticns 

The stability equations for infinitesimal^ two 
dimensional disturbances are formulated under the conventicnal 
assumption that the fluid motion can be decomposed into a mean 
flow (whose stability constitutes the subject of the investi“ 
gation) and a disturbance having pressure ^ and velocity 
components u and v in the streamwise and transverse (or 
radial) direction respectively. Thus in the resultant motion, 
the velocity field is 

u = u + u, v = V + t7, (2.1) 

and the pressure field is 

p=P+p, (2,2) 

where u and V are the mainflow velocity conponents in the 
streamwise and transverse (or radial) direction respectively 
and P is the pressure at any point in the mainflow. The 
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mainflow velocity field is obtained either by the finite 
difference method or by the linearization method both of which 
are described in Appendix A. 

In formulating the disturbance equations, the 
following assumptions are made; 

(i) the fluid is inconpressible and Newtonian with ccnstant 
properties; 

(ii) the mainflow is steady and laminar; 

(iii) the (duct) walls are rigid and impermeable; 

(iv) there are no body forces; and 

(v) the streamwise velocity component U is a slowly 
varying function of the streamwise coordinate x. 


For such a fluid the equations of motion expressed 


in terms of the coordinate system shown in Figure A=>1 ares 



(2.3) 

(2.4) 

(2.5) 


where x is the streamwise coordinate measured from the inlet 
section and y is the transverse (or radial) coordinate 
measured from the channel centre line (or pipe axis), p is 
the density and the kinematic viscosity of the fluid, t is 
the time, and m = 0 for the channel flow while m = 1 for the 
pipe flow. 
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l£ we now substitute eqns. (2.1) and (2.2) in eqns . 
(2.3) through (2.5), make use of equations of motiai for the 
mainflow which are similar to eqns. (2.3) through (2.5) and 
neglect the product terms of disturbance velocity components 
with their spatial derivatives, we obtain the following set of 
linearized hydrodynamic equations for the disturbance velocity 
components and pressure 


+ VL = 0 , 

3 X y 3 y 
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Equations (2.6) through (2.8) can be dimension- 

less by choosing the half width a of the channel (or radius a 
of the pipe) and the average velocity u^ of the mainflow as 
the characteristic length and velocity respectively and 
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Substituting these in eqns. (2,6) through (2.8)/ 

2 

dividing eqns. (2.7) and (2.8) by u^/a and eqn . (2.6) by 

u /a, we get the following set of linear/ partial differential 
0. 

equations for the disturbance properties in nondimensicnal 
form; 


^ u mv a V 
3 X ^ Y a Y 
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( 2 . 10 ) 


( 2 . 11 ) 


(2.12) 


Equations (2.10) through (2.12) are conbined into a 
single equation by differentiating eqn. (2.11) with respect 
to Y and eqn. (2.12) with respect to X/ subtracting the latter 
from the former and using eqn. (2.10). This yields 


3^1 
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u aJ] 

^ 3X 


+ V 
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3Y 
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3X 
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where 


2 - . 21 ) 
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The tern having the coefficient m(m“l) is always 


zero for m = 0 or 1, Henceforth such terms will be deleted 
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Also, the coistant coefficient of m, if any, will henceforth 
utilize the fact that m = O or 1 . Equation (2.13) thus 
reduces to 


3ri ,, aTi ,, 3T1 9 h 9h 

- ^ riJQ ^ 3m _921 jS-i 

^ 3y“ ^ 3X^ 


(2.15) 


It is well known that for the developing flow in a 
pipe or channel, U(X, Y) and V(X, Y) are slowly varying 
functions of X. To express this slow variation we introduce 
another independent variable x^ along X direction such that 


X^^ = e X , (2,16) 

where e is a small dimensionless parameter which characterizes 

the nonparallelism of the flow; e = 0 implies a truly parallel 

flow. The value of e depends on the type of flow and the 

Reynolds number. Bouthier [ss] determined it from the order 

of the terms expressing the boundary layer effect on the 

velocity field and Gaster [63] obtained it in an iteration 

process. For the boundary layer flow on a flat plate, G-aster 

—1 /2 

as well as Bouthier take e = R ^ . smith [68], on the basis 

of comparison of Caster's [63] results with the experimental 

results of Ross et al. [26] and Schubauer and Skramstad [25] , 

—3 /8 

points out that the choice of e = R ^ gives results which 
are in closer agreement with the experimental ones. Unfor- 
tunately no experimental results for the developing flow in 
a channel are available ^d those for the developing flow 
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in a pipe, obtained by sarpkaya [48], are, as Sarpkaya himself 
admits, not free from the effects of non “ ax i symmetric distur- 
bances and of finite level amplitude of the disturbances. 
Therefore, the experimental data of sarpkaya cannot be used 
to estimate the value of e. since the developing flow in these 

ducts, particularly in the near entry region, can be treated 

- 1/2 

as one of boundary layer type, we take e = r for these 
cases . This choice of e is further supported by the express- 
ions given by schlichting [l7] and Goldstein [85] for the 
velocity at the channel centre line and at the pipe axis 
respectively. 

Though e and R are related, we may treat them as 

independent in the expansions that follow. By doing this 

we are, in effect, solving the problem an e -R plane instead 

of on a single curve and thus the real solution is contained 

ti. 

in the family of ficti/ous extensicns over all e and R. 

K 

Now coming back to the eqns. (2.10) through (2.12), 
we note that the coefficients of u, v and their derivatives 
are functions of and Y only. Therefore, the disturbance 
stream function can be taken as 


1{;(X3_, Y, T) = [jZ^q(X^/ Y) + Y) + . . .] 

where 
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^ 'or o 


(2.18) 



20 


wavenumber, and its imaginary part, is the spatial growth 

rate. Using the following transformation relations 
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between (X, T) and (X^ , ©) planes and denoting Y derivative 
with a prime, we get 
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Substituting these transformed relations into eqn. 
(2.15) and equating coefficients of like powers of e we obtain 
order e°s 
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The eigenvalue problems defined by eqns. (2.20) and 
(2.23) or (2.24) or (2.26) are the familiar orr-Sommerf eld 
problems for the parallel flow, in eqn. (2.21)/ which gives 
first order corrections over the parallel flow theory, we note 
that the first two terms on the right hand side represent the 
effects of the streamwise variation of the amplitude of the 
disturbance stream function/ the third accounts for the 
streamwise variation of the wavenumber and of the spatial 
growth rate, and the remaining four terms represent the 
effects of the transverse or radial velocity component of the 
main flow; the last being the only effect considered by 
Shen et al. [49} for the developing flow in a pipe. 

2.2. Boundary Conditions 

Physical considerations such as no slip at the walls 
or boundedness of the eigenfunctions within the flow region 
determine the necessary boundary conditions for eqns. (2.20) 
and (2.21). We now derive them for the channel and pipe flows 
separately. 

2.2.1, parallel Plate channel 

The boundary conditions for eqn. (2.20) at Y = + 1 
are provided by the 'no slip' requirement at the channel 
walls. However, since the velocity profile of the mainflow 
in the channel is symmetric about the centre line, eqn, (2,20) 
with m = 0 admits eigenfunctions which are either symmetric 
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(even mode) or antisymmetric (odd mode). In such caseS/ it 
is ■■ sufficient to consider the stability of the flow confined 
between the centre line and any of the walls, say the, upper 
one (Y = 1). Then the boundary conditions at Y = =1 are 
replaced by the appropriate boundary conditions at Y = 0. 

The boundary conditions, therefore, are 


I M I 



at Y = 0 , 

at Y" = 1 (no slip condition) 


(2.23) 


for symmetric disturbances, and 

n 

^ - 0 = 0 atY=0, 

° (2,24) 

0o = 0 '^ = 0 at Y = 1 

for antisymmetric disturbances. 

The work of Grohne ^94] suggests that for plane 
Poiseuille flew only symmetric disturbances are likely to 
lead to instability. Accordingly, Chen and Sparrow [is] 
studied the temporal stability of only symmetric disturbances 
for the developing flow in a ch^nel. in the present analysis 
both symmetric and antisymmetric disturbances have been 
caisidered for the parallel flow stability analysis, but only 
symmetric disturbances have been considered for evaluating 
the effects of nonparallelism of the flow on the stability 
characteristics. The boundary conditiais for eqn. (2.21) 
which are similar to those for eqn, (2.20) are 



25 


(2.25) 


I n I 

=0 at Y = O , 

f 

and 0^=0^ =0 at Y = X . 

2.2,2. Circular Pipe 

Boundary conditions at the pipe wall are determined 
by the no slip requirement pnd at the pipe axis by the 
requirement that the velocity components and pressure for the 
disturbance be bounded and continuous . Thus we can write the 
following boundary conditions for axisymmetric disturbances 
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boundary condit. 

are 



= 0 

at 

Y = 0 , 

and 

e'l 


= 0 

at 

Y = 1 . 


(2.26) 


( 2 . 21 ) 


(2.27) 



26 


Chapter 3 

SOLUTION OF DISTURBANCE EQUATIONS 

Equation (2.20) together with the appropriate 
boundary conditions depending upon the type of disturbance 
and the flow govern the propagation of an arbitrary, infini- 
tesimal disturbance through an incompressible, viscous, 
parallel flow in a channel or pipe. The variable X- appears 
implicitely in this equation. For given values of w / R 3nd 
U(X^/ Y), this problem can be solved numerically to determine 
the eigenvalue k^(X^) and the eigenfunction P(Y; ) . The 
solution can be expressed as 

jZ^o(Xi, Y) = A(x^) p(Y; X^ ) , (3.1) 

where A(X^) is still an undetermined function. It will be 
determined from eqn. (2.21). If the mainflow is assumed to 
be parallel, A is constant and we have the complete solution 
in eqn .(3.1). 

Substituting for 0 ^ from eqn. (3.1) into eqn. (2.21), 

we get 

= RM(A, p, U, V, k^) , 


where 
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M = [B-,|3 + B.(p 


mfi . -] dA 


Y ■)] ^ + 


L-rfi" _ SLg_i 


2 33^^^ " ^ 


+ [ {B.P + ( 


mS' ' " "' 

Y ^®4> d53f + + ®6^ + ^7^ ^ * 


(3.2) 

The inhomogeneous problem consisting of eqns. (3.2) 
and (2.25) or (2.27) has a solution if / and only if, the 
following solvability condition 


/ M p dY = 0 (3.3) 

o 

'k 

is satisfied. Here P is the eigenfunction corresponding to 
the eigenvalue of the adjoint homogeneous problem 


* „* * 
L (P ) = P 


IV 


2k2p* + k^P* + ^2^ 


3m o* . , 2m, 2 3m,. o* 

-2 P - 

Y Y 


2mk 


* 1 


+ (- 


o 3m, p* T/tti M\rQ* 1 2q* q * \ 

— - - iR L (U}^^ - ^^ ) ( P “ k P - —rr p + — — ) 


Y 


(2u'p*h^p*)K^] 


o 


(3.4) 


With the boundary ccnditions 
p - = p 


= 0 at Y 


O / 


(3.5a) 


= O 


at Y = 1 


for the channel flow (symmetric disturbances) and 

k k 

p = p =0 at Y = 0 and 1 (3,5b) 

for the pipe flow. Substituting for M from eqn. (3.2) into 
ecji. (3.3) and noting that 
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-* a II ^ I VC 

/ ^2 3-^ (P - I P ) P « 


mB, 


it I 

/ + ( 2 B 2 + -p)p + (B 2 + — 


mB^ mB 


2^p*i 8 P 

■)P J 53^ 


dY , 


Obtained by evaluating the integral by parts and using the 
appropriate boundary conditions, we get on simplification 

dA 




(3 ,e 


where ik^ = b^ (X^ )/b^ (X^ ) , 


b^Cxp = 


1 „ mB mB * ■ ™2 

^/C(B2 + ^ + B^)P +( 232 +— )P +^2!^ 


+ / [(B^P - 

o 


4+ - ^ P + B3P ) + B5P + B^p + B^P + vP .]P 


5^(x^) = - / [b^P + B2(P - ^)] p dY . 


The soluticn of eqn. (3.6) is 


A(X^) = A^ exp[i / k^(X^)dX^] = A^ exp[i6/ k^(X^)dx] , 

where A^ is an arbitrary constant of integration. Thus, to 
the first approximation 


,|; (X^,Y,T) = A^P (Y; X^ ) exp [i / (k^ + ek^)dX - iw t] , 


dY , 


(3.7) 


(3.8) 
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where B and k are calculated at each axial location as if 
o 

the basic flow were parallel, and contains the effects of 
the streamwise variation of the mainflow, the eigenfunction 
B, and the eigenvalue k^. To determine , we need to 

evaluate 9 B/3X^ and dk^/dx^ . To do this, we replace by 
P in eqn . (2,20) and the corresponding boundary conditions 
and differentiate the result with respect to X^ to obtain 


R dk 


(3.9a) 


with boundary conditions 



(JA) 

= (■^) 

= 0 

at Y = 0"^ 


(for channel) 

/ 

and 

2Ll_ 

3 X^ 

II 

= 0 

at Y = 1 



(3.9b) 

and 

Li. 

3x^ 

= (Li.) 

0 at 

Y = O and 

1 

( f or pipe ) , 


where 









iRk 

o 

1 

[(V - ^ A,/ 

Y 

3mV 

+ k^v' )p » 

tl 

(P 

-2f)(v' + 

(3.9c) 

and 

^2 = 

4k^(p" - 







iR[u(p" - 2l|-) + B(2(o k^ - 3uk^ =. u" + (3.9d) 


Equations (3.9) have a solution if, and only if, the solva== 
bility condition (eqn. (3. 3)) with M replaced by the terms cn 
the right hand side of eqn. (3.9a) is satisfied. This yields 
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dk 1 ^ 1 ^ 

^ = » [ / P dY]/[/ A 2 P dY] . (3.10) 

o o 

Knowing dk^/dx^ from eqn.. (3.10), 9 P/3X^ can be evaluated 
from the integraticn of eqns. (3.9). 

For a truly parallel basic flow, P is a function of 
Y only, k^^ = 0,. and k^ is a constant. Hence the growth rate 
of any disturbance quantity, such as the velocity, the 
pressure, and the kinetic energy, is given uniquely by the 
imaginary part of k^. On the other hand, the effects of 
ncnparallelism are to make k^ a function of x^, to produce a 
correction ekj_(X^) to k^, and to make the mode shape p vary 
in the streamwise direction. Hence the streamwise variation 
of each flow quantity depends on its distance from the channel 
centre line or from the pipe axis. Moreover, at each distance 
Y, the different flow quantities vary differently in the 
streamwise direction, we shall discuss this point in 
detail in Chapter 4.. 

For determining the effects of nonparallelism of 
the flow on the wavenumber and growth rates, we need to 
integrate the Orr-Sonmerf eld eqn, (2.20) with the boundary 
conditions (2.23) or (2,26). However, for reasons discussed 
in Section 5,1, we solve a set of three coupled, homogeneous 
eqns. (3.12) through -( 3 .15) instead of eqn. (2.20). These 
equations are obtained for an infinitesimal disturbance of 


the form 
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[u/V,p] = [u(Y)/ v(Y)/ p(Y)] exp [i(k^X -uT)] - (3.11) 


by making substitutions in eqns. (2.10) through (2.12). This 
y ie Id s 


- mv . , - 

V + ~ + ik u 
Y o 


= O , 


mu 

u + — rr 


Fu =• RU V - ik^Rp 


= O , 


Rp + Fv + ik u = O , 
o 


2 

where F = k +iR(kU"W). 

o o 

The boundary ccnditions on the basis o£ considerations 
in Section 2.2 are 


(3.12) 


(3.13) 

(3.14) 

(3.15) 


given 


u=p = 0 at Y = 0/ 

u = v = 0 at Y = 1 


(3.16) 


for symmetric disturbances in the case of channel flow^ and 


u = v = 0 at Y = 0^ 

u = v = 0 at Y = 1 


(3.17) 


for antisymmetric disturbances in the case of channel flow 
or axisymmetric disturbances for the pipe flow. Equaticns 
(3.12) through (3.15) may be integrated with proper boundary 
conditions by using any step-by-step integration method; 
the details of the procedure are given in chapter 5. The 
eigenfunction P and its derivatives with Y are then deter- 


mined fron 
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p = iy"' u . p" = u + mu , 

° (3.18) 

p" = [Fu + Rvu + ik^Rp] y^ + mu 

With k , P and its derivatives known, a procedure similar to 
the one above is used to solve the adjoint equations 

(3.19) 

I -k 

ik (v + = 0 , (3.20) 

o Y 

(3,21) 


I * 

* mV • V * r^ 

V + -TT- + Ik U = 0 , 

Y O 

I 

n k 

k rnn 

u + — — “ Fu + ik Rp 


* * ' * 

Rp - Fv + RU U = O , 


with the boundary conditions 

k 

u = p = 0 atY = 0 (for channel flow; 

symmetric disturbances), 

1 

'At 

u = v =0 atY = 0(for pipe flow), (3.22) 

k k 

and u = v =0 aty=l(for both flows). 

k 

The adjoint eigenfunction P and its derivatives with Y are 
obtained from 

I II I 

p = iY^ V /k^, P = Y u , and P = Y u + mu . 

(3.23) 

Since the eigenvalue for the original and adjoint problems 
is the same, no iteration is necessary while solving eqns. 
(3.19) through (3.22). This , in f act , serves as a check 
on the accuracy of the calculated eigenvalues. 
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Chapter 4 

GROWTH RATES AND WAVENUMBERS 


Herein, we derive the relations for the growth rate 
and wavenumber of the disturbance while accounting for all 
nonparallel effects. 

To the first approximation, i.e., to 0(e)/ the 
stream function of the disturbance is given by eqn. (3.8), 
which may be rewritten as 

‘4' = A^P (Y; X^) exp [i / a dx = icu'T] , (4.1) 


where a * k + e k^ . 

o i 

The amplitude of this stream function is 
!4)| = I ^ ^ ^1^ I exp(- / dX) , (4.2) 


where is the imaginary part of a. We define a growth rate 
based on '|> in x-space as 


= 1 aUn 


(4.3) 


Using eqn. (4.2), eqn. (4.3) becones 






3 I PI 
a x^ * 


(4.4) 


we note that the growth rate of the stream function is 
dependent on the streamwise as well as the transverse (or 
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radial) coordinate since P is a function of both and Y. 

Alsor when k . = O, i.e.y at the neutral points determined 
oi 

by the parallel flow theory, there is still growth or decay 
due to the higher-order effects. Thus, the higher-order 
corrections are essential in determining the correct neutral 
points for the stability analysis of nonparallel flows. 

Another surprising feature of 0(e) corrections is 
that the growth rate is different for different flow quanti- 
ties. Consider, for example, the x-component of the 
disturbance velocity, given by 


u 


1 _ ^ 
^ a Y 


^ 9 R 

exp(i / a dx - iw T) . 


(4.5) 


The growth rate of u is then given by 


u 


a. 4 - 



3 1 8 p/a Y 
a X, 


(4.6) 


This, in general, is different frcsm the growth rate of the 
stream function, given by eqn. (4.4). The question then 
arises about the "proper measure" for the growth of the 
disturbance, if one has experimental data to compare with 
the calculations, one must obviously use the same quantity 
as that which was observed. Unfortunately no experimental 
results for the stability of developing flow in a channel 
are available and those for the pipe flow, obtained by 
Sarpkaya [48] , are, as saiqskaya himself admits, not free 
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from the effects of non^axisymmetric disturbances and of 
finite level amplitude of the disturbances. Moreover^ 
though sarpkaya mentions that the growth rate of the stream^ 
wise component of the disturbance velocity was measured at 
different radii, he neither reported their magnitude nor the 
radii at which measurements were made, since these details 
could not be obtained [95] , it was decided to compute all 
the growth rates mentioned above as well as the growth rate g^ 
of the kinetic energy density of the disturbance, defined as 


gg(X) 


.1 ^ 

dX 


«i 


£ a c 
2c 3X ' 


where e 


I (2ttYf (u^ + v^) dY , 
o 


and 



^ + I kj ^ I Pi dY . 


(4.7) 

(4.8) 

(4.9) 


Here E represents the mean kinetic energy density, averaged 
over time and integrated across the flow domain . For the 
growth rate based on E, we have included a factor of half in 
the definition so that comparison with other growth rates is 
possible . 

From eqn . (4.1) one notes that not only the growth 
rate but also the wavenumber is affected by the ncnparallelism 
of the flow. The modified wavenumber for any disturbance 
property G is obtained from 


Wavenumber 


3 [arg (G)] /ax 


/ 


(4.10) 



36 


where arg(G) represents the phase of the complex quantity G. 
One, therefore, gets different wavenumbers for different 
disturbance properties. The modified wavenumbers ^#1 (X, Y) 

and Y) based on the stream function and streamwise 

component of velocity of the disturbance are respectively 
given by 


(X, Y) = + lx / (4.11) 

3 

VJN^(X, Y) = +3“x , (4.12) 

where k^^ is the real part of k^ . It was found that the 
effect of nonparallelism of the flow on the parallel flow 
wavenumber, insignificant since k^^ is much larger 

than ek^^ and the terms representing the effect of streamwise 
variation of the argument of the eigenfunction. Therefore, 
the modified wavenumbers, though computed, are not reported 


here 
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Chapter 5 

COMPUTATIONAL PROCEDURE AIsfD DETAILS 

In this chapter we first describe a compariscn of 
the various methods available for controlling the parasitic 
errors accruing during numerical integration. This is followed 
by details of the method found to be most efficient and 
ecoionical and a brief outline of a technique to search for 
eigenvalues on the complex k^=plane. The chapter is closed 
with computational details. 

5.1 . conparlson of various Methods to Control the 
Parasitic Error 

Equations (3.12) through (3.14) with (3.16) or (3.17) 
can be integrated by any of the suitable integration schemes 
available, in the case of pxpe flow, L' Hospital rule is ueed 
to simplify the terms involving Y in the denominator for 
starting the integration of these equations. One can take 
advantage of the linearity of the equations to change the 
two point boundary value problem to an initial value problem 
using complementary functions [96]. Since two boundary 
conditicns are provided at Y = O as well as at Y = 1, we get 
two solutioi vectors irrespective of the end from which 
integration may start. If we start integration from Y = O for 
given values of R, w and an initial guess for with the 
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initial linearly independent and orthonormal solution 
vectors, (1, 0, O, O) and (O, 0, 1, 0) for the channel or 
(O, 1/ O, 0) and (O, O, O, 1) for the pipe flow stability 

r=, — 

corresponding to (v, u, u , p), we get two solutions which 
may be represented by the subscripts 1 and 2. The boundary 
conditions on the wall demand that the determinant 


D 


u^(l) 

v^(l) ^2^^^ 


( 5 . 1 ) 


must vanish. This condition determines the eigenvalue k^. 

The above method, though mathematically exact, 
fails to give any meaningful results when applied as a nume= 
rical procedure due to the parasitic growth of one of the 
two independent solutions, specially when the Reynolds number 
is large, say a few thousand, as is true in most of the 
hydrodynamic stability problems including the present ones. 
For this reason Kaplan [97] filtered out a part of the 
growing solution from the slowly growing solution at each 
step of the integration. Davey [?] cleverly circumvented 
the situation by dividing the integration range into a nurnheir 
of steps, say c, integrating the given differential equation 
separately for the range of each step starting everytime 
from the same set of linearly independent and orthcnormal 
vectors and finally combining the solution vectors obtained 
at the end of each step which he called as transfer 
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matrices. Godunov [98] and Bellman and Kalaba [99] indepen™ 
dently used the Gram -Schmidt orthcnormalization procedure to 
make the solution vectors linearly independent during 
numerical integration. This method has also been used and 
discussed by ccnte [lOo] , wazzan, Okamura and smith [ss]/ 
Roberts and Shipman [96], Scott and watts [lOl] , and 
Huang ri02] , 

we tried all the above techniques with the fourth 
order Runge-Kutta integration scheme for the channel flow 
since this was the first flow geometry considered during 
computation. We also used the method described by Lee and 
Reynolds [l03] . Table 5.1 shows the approximate time 


Table 5.1. Conparison of different methods. 


s. 

M o « 

Method of 

1 

JTime/ 
[iteration 
[(second ) 

Upper limit of Rk 

t'.t' Qj- 

1. 

Lee and Reynolds [l03] 

30 

Not determined 

2. 

Davey [?] 

5.5 

60000 with c = 1 

3. 

Kaplan [97] 

5.4 

28000 

4. 

Gram-Schmidt [lot/1 
for eqns. (2.20) with 
(2.23) 

3. 3-3. 5* 

Convergence becomes 
slow at high R .and oj 

5. 

Gram-Schmidt for 
eqns. (3.12) to (3.14) 
with (3.16) 

2. 5-3. 2* 

No limit 


* Actual value depends on the number of orthcnormalizations. 
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required by the different methods for one iteration on IBM 
7044 conputer in double precision mode with the step size 
of 0.01; eqns. (2.20) with (2.23) being integrated. Also 
given in this table is the time required per iteration while 
integrating eqns. (3.12) through (3.14) with (3.16) as a 
set of three simultaneous equations using the Gram-Schmidt 
orthonormalization procedure as well as the approximate upper 
limit of the product upto which the different methods 

worked satisfactorily. 

Lee and Reynolds method was used in the initial 
stage and was soon given up as it takes sufficiently large 
computation time and convergence to the eigenvalue is very 
poor. Davey proposes that his method can be used upto 
very high values of Reynolds number by choosing c properly. 
However, in comparison to the Gram-Schmidt orthonormalizaticn 
procedure, it not only takes more time but also leads to a 
highly involved procedure for determining the eigenfunctions. 
Kaplan's filter technique besides requiring more computation 
time suffers from two disadvantages. These are: 

(i) It fails when the range of integration contains more 

than one region in which viscous forces are important, 
as pointed out by Sharma [los] . Thus in the case of 
pipe flow where there are two viscous regions, one 
near the boundary wall and cne away from the wall/ 
the integration will be inaccurate over the wide range 
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between the two viscous regions, sharma suggests that 
whatever filter one uses in Kaplan's technique in such 
cases, integration through a viscous region will be 
inaccurate. Davey and Nguyen [106] found that for 
the central mode in pipe flow stability analysis, 
Kaplan's method did not work at all. 

(ii) It cannot be used at high Reynolds numbers and it is 
rather cumbersome to use for a higher order differ- 
ential system when a multiple filtration is required. 

From Table 5.1 it is observed that the Gram-Schmidt 
orthcoormalization procedure is the best possible method 
ammgst the various methods available. This view is also 
supported by Gersting, Jr. and Jankowski [l07] , which came 
to our notice after our tests with various methods were over. 

Table 5.1 also reveals that integration of eqns. 
(3.12) through (3.14) instead of a single Orr-Sommerf eld 
eqn. (2.20) is economical, it was observed that for the 
given values of u , R and the same initial guess for k^ 
there was a significant decrease in the number of points at 
which orthonormalization was necessary as we switched over 
from the integration of Orr-Sommerf eld equation to the 
integration of the set of eqns. (3.12) through (3,14). This 
implies less number of ccmputations, less round off error 
and hence better accuracy. It was also observed that for 
the same starting value, the number of iterations required 
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for convergence to the true eigenvalue is decreased, thus 
requiring less time to compute the eigenvalue when eqns. 

(3.12) through (3.14) are integrated. In a recent publi- 
cation, Antia [lOS] points out that in hydrodynamic stability 
problems the elimination of variables to get higher order 
equations introduces singularities in the resulting equations 
and so the numerical treatment becomes difficult. 

The results reported here for the channel as well 
as pipe flows were, therefore, obtained by integrating eqns. 
(3.12) through (3.14) using the Gram-Schmidt orthonormalization 
procedure described briefly in the next section along with 
the improvement for iteration to the eigenvalue. This 
enabled us to determine the eigenvalue f or Reynolds number 
upto 10 without any difficulty. 

5.2. Gram-Schmidt Orthonormalization Procedure 

The basic idea in this procedure is that during the 
integration the base vectors are orthonormalized each time 
the vectors start to lose their numerical independence as 
judged by a specified orthonormalizaticn criterion. The 
orthcxiormalized solutions, defined on various subintervals , 
are then pieced together to obtain the desired soluticn. 

In the following sections, we discuss the procedure, deter- 
mination of eigenfunctions, orthonormalization criterion 
and the application of the method to an eigenvalue problem. 
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5.2.1. Procedure 

we choose n steps of uniform mesh size Ay for the 

region 0 ^ Y ^ 1 and apply any standard integration method 

1 2 

to obtain the base solutions S(Y) = [s , s ] at the mesh 

points Y / Y-, ; the superscript denoting the column 
. o 1 n 

1 2 

number of the matrix. For example, s and s may be taken 
as 

s^ = [v.y u., u., , (j =1,2) . (5.2) 

~ J J J J 

we examine these solutions at each mesh point as we integrate 
and when the base solutions exceed certain orthonormalization 
criterion (cf. section 5.2.3) we orthonormalize the solutions. 
Let Y^ be any such point at which orthonormalization is 
first done. The matrix of solutions S(Y. ) is orthonormal- 
ized to get Z(Y.) where 2 and s are related by 

z(y^) = S(Y^) ; (5.3) 

being a nonsingular upper triangular matrix with elements 

b. .. Given the matrix S(Y.), we can obtain the matrices Z 

13 *» 1 ~ 

and B by applying the Gram=Schmidt recursion formulae 

12 

(eqns. (5.4)) for orthonormalizing a set of vectors [s , s ]. 

mm 

Using ( , ) as the usual notaticn for inner product of 

two vectors, we get 

b - (s^ sh“^/2 - b s^ 

°ii ‘‘f ' ® ^ ' 5 ", hi - ' 
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2 ,2 1 , 1 

w = s •=(s,z)z, b 


22 


" ^ 22 ^! ' 


(5.4) 


b 


12 


< 5 ^' bl ’"22 ' *=21 ' ° • 


It is clear that the elements b^^ are available as 
a byproduct during the orthonormalization process. iSTow we 
take Z and not S as the starting solution for the next step 
and continue integration. 

5.2.2. Eigenfunctions 

Due to orthonormalization during integration / the 

base solutions are changed several times. Let us use the 

1 2 

notation Q(Y) = [q / q ] to denote the' discrete solution 

available as a result of this process. Let us also suppose 

that we encounter b points at which orthonormalizaticai is 

carried out and that Y = 1 is the bth orthonormalization 

n 

point. Note that b n. At Y = Y^ / orthonorm alizaticai gives 

Z(Y„) = Q(Y^) . (5.5) 

^ n n ^ 

The function S(y ) at Y = Y is now obtained from 
- n n 

where is the solution of the system of equations obtained 
from the boundary conditions at Y = 1. At other mesh points 
the function is obtained by backward transformation in the 
following manner. The matrices B at all the b 
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orthonormalization points are stored during the forward 
integration. Then the funcuicn S(Y^) ^t any point i between 
the two consecutive orthonormalization points j and (j-1) 
is given by 

S(Y. ) = Q(Y )y^ , 1 ^ ' (5.7) 

where will be a constant vector for all y^ between the 
two consecutive points of orthonormalization. As soon as 
the point i moves to the left of the orthonormalization 
point (j-1)/ the vector y ^ recomputed from 

_^(j“l) = j = t)/ b-1 / .../ 1 . (5,8) 

The backward resolution is continued until we arrive at the 
initial point Y = 0, where an immediate check on the accuracy 
of the resolution is available by comparison with the boundary 
conditions at Y = O. 

5.2.3. Orthonormalization criterion 

In the method discussed above/ the orthonormal- 
ization at each mesh point is neither necessary nor desirable 
due to loss of some accuracy by too frequent an orthonormal- 
ization. several criteria for determining when orthonormaliza- 
tian is to be performed are available [lOO, lOl], The 
easiest criterion to implement and very inexpensive to 
compute is the one proposed by Conte [lOO] that requires an 
orthonormalization to be carried out whenever the magnitude 
of any vector S(y) exceeds a preassigned constant K, that 
is, whenever 
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max (s^, >K/ (j=l,2). (5.9) 

j ~ ~ 

The constant K depends in some way upcn the magnitude o£ the 
expected solution and is, therefore, highly problem- 
dependent as pointed out by Scott and >Jatts [lOl]. However, 
for the present problems, we found that taking K as 100 
works very well. 

5.2.4. Application to Eigenvalue Problems 

While discussing an eigenvalue problem 'involving 
the Schrodinger equation, Conte [lOo] suggested that the 
orthonormalization be the same for all iterations to the 
eigenvalue. He, therefore, proposed that the set of b 
orthonormalization points and the corresponding matrices B 
should be determined for the first choice of the eigenvalue 
and thereafter the same matrices at the corresponding points 
should be applied for all successive iterations, while this 
is a sound suggestion in order for the successive approxi- 
mations to the eigenvalue to be consistent, it works only 
if the eigenvalues of the system of differential equaticns 
are not too widely separated. For hydrodynamic stability 
problems, it implies that Conte's above suggesticn fails 
for relatively high values of the Reynolds number. What 
happens is that if the matrices B (even at the same set of 
b orthcnorraalization points as for the first choice of the 
eigenvalue )• are determined for successive iteraticns. 
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they are found to change from one iteration to the other. 

Thus/ if the same matrices B are kept for successive iter ac- 
tions as for the first one, proper orthonormalization of the 
solution during successive iterations is not possible 
resulting in "overflows". However, if the matrices B are 
determined anew for every iteration, proper account for it 
will have to be kept so that successive approximations to 
the eigenvalue are consistent. The following method is 
therefore adopted. 

For an eigenvalue problem, the starting value, of 
S(0) is taken as 

(5.10) 

where let be a (2 x 2) null matrix and S 2 be a (2 x 2) 
identity matrix. Every orthcnormalization changes the 
initial matrix S(0) by making an upper triangular matrix 
[ 99 ], while remains a null matrix, we are concerned only 
with the changes in the values of the diagonal elements Sjj 
of s.p. clearly, s.. is, in effect, multiplied by b ., 

(j = 1, 2), each time an orthonormalization is carried out. 

The boundary conditions at Y = 1 are generally satisfied, 
for an eigenvalue, by setting an appropriate determinant, 
such as D in eqn. (5.1), to zero. The value of this deter- 
minant depends on the initial values S(0) and thus cm the 

values of the diagonal elements s.. of S„- Therefore, for 

33 ~ z 
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the successive iterations to the eigenvalue to be coisistent, 

the determinant D should take into account the changes in 

s.. due to orthonormalization. This is achieved by dividing 
33 

D by the product o£ all bj^z (j =1/ 2), for all the b 
orthonormalizations during any iteration. Thus if 

b 2 

' N = n n b":. z (5.11) 

i=:i j=i 

where the superscript on b^^ represents the latter's value 
for the ith orthonormalization z the determinant D/Kf will 
correspond, during each iteration, toSjj =1, (j =lz 2), 
and successive iterations will be consistent. 

while the above technique is mathematically 
correct, a practical difficulty may arise in canputing the 
value of N since sane, if not all, of the b values of b^^ 
are of order (K ^). A direct use of eqn. (5.11) to find N 
may, therefore,, result in "underflow" on the computer. To 
avoid this, N was found from 

N = n W I b^ . z (5.12) 

i=l ° 

where N is a fixed constant of -order (K). Like K, N is 
o o 

dependent cn the problem. For K = 100, = 140 for the 

channel flow and = 100 for the pipe flow worked very 
well at all axial locations, Reynolds numbers and frequencies 
except at X = x/R = 0.0005 for the pipe flow where was 


taken as 130. 
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vihile using eqn. (5.12), we must make sure that the 
number of orthonormalizations, b, is the some for every 
successive iteration to the eigenvalue, or that is used 
the same number of times during every iteration. If it is 
not so, proper account must be kept in order to make succe- 
ssive iterations meaningful. In fact, the set of b ortho- 
normalization points may be found for the first iteration 
and thereafter kept same for the successive iterations. 

5.3. Eigenvalue search Technique 

For the pipe flow case, a few eigenvalues of the 
Hombeck profile for different R and w at a given axial 
location were determined by using the eigenvalue search 
technique [40] . Although a ccmplicated technique for 
isolating a number of these eigenvalues simultaneously 
rather than one at a time b^s been developed by Delves and 
Lyness [l09] , it was not used for the present problem 
because the main purpose here was not to find a number of 
stable eigenvalues but to determine the unstable eigenvalue, 
if any. 

The search technique follows directly frcm the 
argument principle in complex -variable theory. For a 
complex function f (z) , ]| analytical except for poles in the 

interior of a closed ci|rve c in the complex z— plane and for 

|| 

f(z) and its first derivative continuous cn C, Cauchy's 
theorem gives 
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N - P 


2% i 


a; 

f (z) 


dz j. 


or N-P = 2'n;i m/2h; i = m / 


(5.13) 


where N and P denote respectively the number of zeros- and 
poles of f(z) within the closed region C (counted with their 
multiplicities) and M is an integer denoting the net multiple 
of 21^ by which the phase angle of f(z) changes as z moves 
once around on O . 

For either flow stability, each element of the 

determinant D in eqn. (5.1) is a function of the complex 

eigenvalue k^, and this function has no poles [40], that is, 

P = O in eqn. (5.13) and N = M. Thus the problem of 

finding the number of eigenvalues within a closed region of 

the k^-=plane is equivalent to counting the net multiples of 

211 by which the phase angle of the determinant changes as 

k assumes values cn a closed contour in the k -plane, 
o o 

For -the stability of the Hornbeck profile in the 
pipe, sane regiois of the fourth quadrant on the k^-plane 
were examined by the above technique for possible (unstable) 
eigenvalues at a few selected values of u , R and X. If 
an eigenvalue was indicated, it was bracketed by dividing 
the region into smaller subregions until Muller's iter- 
ation converged to the true eigenvalue. For other stability 
problems, use of this technique was not required as some 
starting values were available. 
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5.4. ccxnputational Details 

In order to select an appropriate step size for the 
fourth order Runge-Kutta method some eigenvalues were 

determined for the channel flow at different X, and at each 
X for a few combinations of R and w with three step sizes 
of AY = 0.01/ 0.005 and 0.0025. The error in the eigen- 


value with the smaller step size was then calculated from 

[no] 


error in eigenvalue 
with 

the smaller step size 



(5.14 ) 


where L is the order of the integration method and 

and (k ) are the eigenvalues determined with the step 

sizes ^Y and AY/2y respectively. The error in the eigen- 

-5 

value with AY = 0.0025 was found to be of the order of 10 
or less in both the real and imaginary parts of k^. This 
error was much less than that in the eigenvalue obtained 
with AY = 0.005. It was also observed that the computation 
time per iteratiai was almost doubled as the step size was 
halved. Therefore^ a compromise between computer time and 
accuracy was made and a step size of 0.0025 was selected. 
Muller's method [ill] was used for iteraticn to the 
eigenvalue. It is an iterative procedure for finding the 
real and complex roots of a polynomial equation whose 
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approximation to the root is taken to be one of the zeros 
of the second degree polynomial which passes through the 
functional values corresponding to three selected points. 

|7 

The iterative procedure is continued by droj^ing the first 
point and considering the second point as the first point, 
the third point as the second point and the new point 
obtained above as the third point. The advantages of this 
method ares (i) the canplex roots can be obtained; (ii) the 
iteration requires only the evaluation of functional values 
and does not include the value of the derivatives of the 
function; and (iii) after the three initial estimates have 
been processed, only a single pass through the integration 
procedure is required for each iteration. 

Since the introduction of the factor N for 

o 

maintaining consistency between successive iterations makes 
the determinant D (eqn. (5.1)) dependent on and hence 
arbitrary to seme extent, the following convergence criterion, 
instead of making D less than a preassigned value, was used: 



where Ak is the difference between the eiaenvalues of the 
o 

last two iterations and 6 is seme preassigned censtant* 

For the accuracy mentioned above, 6 was taken as lo” . 

In almost all the cases, convergence was achieved in 3-4 
iterations, while determining the neutral points, further 
iteration was terminated when Ik .1 became less than lo"*^. 

* m ' 
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While finding the eigenvalues for the Hombeck 
profile / the eigenvalue search technique described in 
section 5.3 had to be used in a few initial cases since it 
was not possible to achieve convergence while starting with 
the values for the Sparrow profile. 

For calculating velocity profiles in the entrance 

region of the channel and pipe by the Sparrow method/ eqns. 

(A. 14) and (A. 16) respectively were used. The Bessel 

functions and for real arguments, appearing in eqns. 

(Ails), were computed using Garg's programme [ll2] which was 

simplified to suit our requirements. The infinite series in 

eqns. (A* 14) or (A«16) was terminated when the ratio of the 

last term retained to the partial sum of all previous terms 

—18 

became less than lo" . The velocity gradients -were deter= 
mined by evaluating the relations obtained by differentiating 
the eqns, (A.14) or (A, 16). Fifty values of X were obtained 
for both the channel and pipe flow by using Muller's root 
finder [ill] for solution of the characteristic eqns. (A.,15) 
and (A«17). These roots are tabulated in Table Al. , 

For calculating the velocity field at different 
cross-sectiorp by the finite difference method, eqns, (A*-6) 
through (A. 8) were solved with AY = 0.0025 at all cross- 
secticns and Ax = lo”^, lo“^, 5 x lo"*^ for 0 < X < 0.01, 

0.01 < X 0.1 and 0,1 < X < 0.2, respectively. This 
axial spacing was decided cn the basis of our efforts to 
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reproduce the velocity profile given by Hornbeck [lis] for 
both the flow problems. The velocity gradients were deter- 
mined by use of central differences involving an error of 
2 

Oi ). For use during the Runge-Kutta integration the 
velocity at the mid-step points ^ such as the point Q in 
Figure 5.1 ^ was determined by taking the average of the 
velocities at the point Q obtained by quadratic interpolation 
first fron the known velocities at step points A, B and C 
and then from the known velocities at points B, c and D. 

This procedure was 


k AY Ay 

^ — o-H 

! 2 ] 2 ! 

■e e 5€ © — 0 

APBQCRD 1 

Y 

Fig. 5.1 Seme step and mid -step points in the 
range of integration. 

adopted after careful numerical experimentation in which 
velocity at step points like B, D, etc. was calculated while 
treating them as mid-step points. It was found that the 
above procedure resulted in values which differed from the 
actual velocities at these points in the 7th or 8th decimal 
place. Linear or simple quadratic interpolation resulted 
in much larger error. However , at the mid-step point 
(1 - i^Y/2) f the velocity was determined by simple quadratic 
interpolatioi only. 


0 ^ 

0 
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Once k / P/ and 6 were found , dk /dX-, was calcu- 
o o 1 

lated from eqn. (3.10) using fifth order composite Newton- 
Cotes quadrature formula [114] for finding the integrals 
numerically. The values of 3 P/9X^ are then given by the 
integration of eqns. (3,9). it may be noted that for the 
pipe flow case, the integrands in eqn. (3.6) vanish at Y = 0. 
Also L 'Hospital rule is used to simplify the expressiais for 
and A 2 at Y = O in eqns. (3.9c, d) for starting the 
integration of eqn, (3.10). 

For given co and R, eqns. (3.12) through (3.14) 
with (3.16) for channel and (3.17) for pipe case were also 
solved to yield k (X- ) , P and its Y derivatives at three 

O 1 

axial locations x - Ax, X and x + AX; AX being taken as 10 

Using central differences involving an error of O (R AX) , 

6k^/6x^ and 6p/6x^ were then obtained. It was found that 

the values of 6k /6X- and dk /dX, were in agreement within 

o 1 o 1 

computational accuracy. Also dp/Sx^^ ^nd 9p/9X^ were in 
agreement at every point in the flow domain. 

calculations were performed partly on an IBM 7044 
computer and partly on a DEC 1090 computer for the channel 
flow stability and completely on the DEC 1090 computer for 
the pipe flow stability problem. The computation was done 
in double precision mode. Both the computers carry 17 
digits in this mode, A FORTRAN 10 versiai of the computer 
programme for nonparallel pipe flow stability is listed 
in Append iXj D. 
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Chapter 6 

RESULTS AND DISCUSSION 

The computation for determining the spatial 
stability behaviour of the developing flow in a two dimen=- 
sional parallel plate channel and in a rigid circular pipe 
was done in three stages, viz., (i) the determinaticxi of 
velocity field at 401 equi-spaced points along the Y 
direction by the finite difference and linearization methods, 
(ii) the integration of the continuity and mcsmentura equaticns 
for the disturbances superimposed on the developing flow, 
assumed parallel for this purpose, and (iii) the determin- 
ation of various growth rates based cn different distur- 
bance properties for the actual ncnparallel flow. The 
following sections discuss the results in this order. 

Since the velocity profiles were computed at seme preselected 
values of x instead of X, the growth rates, reported herein, 
were also obtained at various X, and hence they are, here- 
after, referred to as functions of X instead of X; X being 
X/R . 

6.1 . parallel Plate channel 

Spatial stability results have been obtained for 
the B-o profile at several values of X and also for the 
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Sparrow profile at three values of X for checking our 
computer programme and for determining the effect of different 
velocity profiles on the stability characteristics of the flow. 

Developing velocity Field 

Figure 6.1 shows the velocity and its first deri- 
vative with Y obtained by the methods suggested by Bodoia 

« 'k 

and osterle and by Sparrow et al, at X = 0.00208 (X = 0.005) 

... 'k 

and at X = 0,00408 (X = 0.009). It is observed that the two 
velocity profiles differ more frcm each other near the 
boundary layer edge but the maximum difference is within 5% 
at X = 0,00208. The velocity gradients also differ consi- 
derably near the channel wall and the boundary layer edge; 
the difference near the boundary layer edge being as much as 
100%. The velocity gradient for the B-o profile is larger 
than that for the Sparrow profile near the channel wall but 
f alls below 'it mid-way through the boundary layer thickness. 
These differences decrease with increasing X. At X = 

k 

0.08393 (X = 0.10) the two methods give almost the same 
velocity field. The B-0 and Sparrow velocity profiles at 
same values of X are shown in Figure 6.2/ and are also 
tabulated in Tables A3 and A4 . 

6.1.2. parallel Flow Instability to Symmetric and 
Antisymmetric Disturbances 

Grohne [94] has shown that for fully developed 
flews having symmetric velocity distribution, symmetric 
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disturbances are more unstable. However, it was decided 
to verify this statement in the case of developing flow 
through a channe 1 . The ref ore , eqns . (3.12) through (3.14) 
with the boundary conditicais (3.16) for symmetric distur- 
bance or (3.17) for antisymmetric disturbance were 
integrated to get eigenvalues for the unstable m.odes for 
different combinations of to and R at several axial loca- 
tions. Figures 6.3 and 6.4 show the neutral curves ( (^ vs. 

R) and(k:^^ vs. R) at X = 0.004, 0.006 and 0.008 for symmetric 
and antisymmetric disturbances. Tables Bl and B2 give the 
values of frequency and wavenumber corresponding to the 
neutral points at different X and at each x for different 
Reynolds number for both types of disturbances. It is 
observed fran these tables that at X = 0.001 and 0.002 the 
difference in the neutral point values of w and at any 

R for symmetric and antisymmetric disturbances is too small 
to be displayed graphically. As X increases the difference 
in the two neutral curves increases and their lower 
branches cross each other meaning thereby that there exist 
seme eigenstates which are more unstable for antisymmetric 
disturbances. However, the critical Reynolds number is 
smaller at all x for the symmetric disturbances. Further, 
the difference between the critical Reynolds numbers for 
symmetric and antisymmetric disturbances increases as X 
increases. This is quite evident from Figure 6.5 which 
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shows the variation o£ R / w and k with X for both types 

c c or^ 

of disturbances, it is clear that the difference between 
critical frequencies and wavenumbers for the two types of 
disturbances also increases with X. The critical Reynolds 
number for symmetric disturbances continues to decrease while 
that for antisymmetric disturbances passes through a minimum 
value of about 12100 at X = 0.005 and then increases rapidly. 
It iS/ therefore^ expected that the fully developed flow 
will be stable to antisymmetric disturbances. It is thus 
obvious that the developing flow in a channel is more 
unstable to symmetric disturbances in comparison to anti- 
symmetric disturbances except for a few eigenstates in the 
near entry region. However/ we need consider only the 
symmetric disturbances so far as our interest is in finding 
the critical values. 

6.1.3. Effect of velocity Plstribution on Parallel Flow 
Instability — Symmetric Disturbance's 

The variation of critical Reynolds number with X 
and with x for the present analysis of the B-0 profile and 
for Chen's temporal stability analysis [lis] of the Sparrow 
profile using the finite difference scheme of Thomas [s] 
is shown in Figures 6.6 and 6.7. it is seen that the 
critical Reynolds number for the two profiles differ 
significantly in the region close to the entrance plane 
where the B-0 profile gives a lower critical Reynolds number. 
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At X = 60 the critical Reynolds number for the B-O profile 
is about 10900 which is almost half of that for the sparrow 
profile. This difference between the values decreases as 
X (or X) increases to X - 0.084 (or X =; 440) where the two 
curves appear to coincide with each other and remain so 
thereafter. Such a behaviour is to be expected since the 
two velocity profiles, as noted above, merge into cne at 
X = 0.084, One, therefore, draws the conclusion that the 
larger difference in critical Reynolds number in the near 
entry regicxi is due to the difference in the two velocity 
profiles. The three points marked ® on Figures 6,6 and 6.7, 
obtained for the sparrow profile by the present technique, 
serve as a very good check for the present method, Chen 
fll5] used several methods including the finite difference 
method of Thomas [3 ]. Though Chen himself regarded the 
finite difference method to be most accurate, it is strange 
that very few results reported by him were obtained by the 
finite difference method. 

Figure 6.8, where the neutral curves for the 
present analysis at two values of X (both for the B-0 and 
Sparrow profiles) as well as those of Chen obtained from 
the regular viscous solutiois for the full channel profile, 
i.e., the Sparrow profile, have been brought together, 
shows further the effects of the difference in the two 
velocity profiles on the stability characteristics of the 



61 


flow. The data for the Sparrow profile are given in Table B3. 
The neutral curves for the B-0 profile encompass more 
unstable region. Therefore, the flow with the B-o profile 
is unstable for a wider range of frequencies at a given 
axial location and Reynolds number. Figure 6.8 also shows 
the difference in the neutral curves for the same (Sparrow) 
profile obtained by two different methods. Recalling the 
above discussion on Figures 6.6 and 6,7, we note that our 
results, even for the Sparrow profile, are more accurate 
than those of Chen shown in Figure 6.8. Chen [HS] gives 
several tables of data for the neutral curves obtained by 
the asymptotic method of Heisenberg [ll£] , viz,, the regular 
viscous and composite solutions for the full channel 
profile, and for the velocity profile obtained by treating 
the developing flow as a boundary layer. All of these 
neutral points are less accurate than those obtained by the 
present analysis or by the finite difference method used by 
Chen. Chen's neutral curves obtained from the regular 
viscous solution for the full channel profile are found to 
be closest to the present results, and are, therefore, 
included in Figure 6.8 for conparison. Besides, Chen's 
analysis is for temporal stability while the present work 
considers the more realistic spatial stability of the flow. 
Figure 6.9 shows the neutral curves vs, r) at X = 

0.00208 and 0.00408 for the B-0 and Sparrow profiles 
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obtained by the present analysis. The same trend as noted 
for Figure 6.8 is observed. 

The neutral curves/ ( w vs. R) and vs. R)/ at 

various values of X are shown in Figures 6.10 and 6,11. It 
is observed that the Reynolds number, frequency and wave- 
number at the critical point and area of the unstable region 
decrease with increasing X. The slope of the lower branch 
of the neutral stability curve is very large at X = 0.001 
and it decreases with increasing X whereas the upper branch 
of the neutral curve is almost flat at all X. The critical 
Reynolds number for the fully developed flow was found to 
be 3848.1/ which agrees very well with the value of 5772.22 
reported by Davey [ll7] and orszag [s] , after multiplying 
it by the factor of 1.5, since we have used average velocity 
in place of the maximum velocity as the reference velocity 
for nondimensionalization . variaticn of the critical wave- 
number and frequency with X is shown in Figure 6.12. One 
notes that the critical frequency and wavenumber decrease 
with increasing axial distance in the entrance regicn and 
approach asymptotically the corresponding values for the 
fully developed flow; the two curves being nearly parallel. 

Figure 6.13 shows a representative variatioi of 
kor vs. to at X = 0.001, 0,002, 0.004 and 0.008 for R = 

16000. and at X = 0.001 for R = 18000 and 20000. One notes 
that the different curves are almost parallel. This 
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behaviour helps a great deal in determining the guess value 

of so that ccnvergence is easily achieved. 

Figure 6.14 displays the real and imaginary parts 

of the eigenfunctions for the streamwise and transverse ccm'= 

pcnents of the disturbance velocity at X = 0.001. These 

have been normalized with respect to their respective maximum 

absolute values; the ratio of iOl to 1 vi being 13,93, 

' 'max 'max 

It is obvious that the boundary conditions (3.16) are 
satisfied . 


6.1.4. Effects of Ncnparallelism of the Flow 


Growth rates based an ip and E were obtained at 
X = 0.001, 0.002, 0.004, 0.006 and 0.008. Figure 5.15 shows 
g^ , the growth rate based on ^ , as a function of Y f or 
different combinations of X, R and u. It is clear that in 
the region near the channel wall the dependence of on Y 
is quite strong. It is also observed that the maximum growth 
rate occurs at the centre line of the channel, while g^ 
decreases gradually and uniformly with Y upto the boundary 
layer edge for all combinations of X, R and o) , its varia- 


tion beycnd this point depends on the values of X and upon 
the position of the (o) , r) point relative to the neutral 
cuirve on the w - R plot. If the selected combination of R 
and (0 lies close to the neutral curve, g^ decreases suddenly 
near the boundary layer edge and then increases near the 


channel wall; the magnitude of the depression decreases as 
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X increases or as one goes into the stable region away from 
the neutral curve. Maximum depression in the growth rate 
curve at any X and R is found to occur at frequencies midway 
between those corresponding to the upper and lower branches 
of the neutral curve. The location of this depressich shifts 
with R and w at a fixed X. 

Since g^jj is maximum at the centre line, it, g^ (X,0) 
is compared with other growth rates such as that obtained 
from the parallel' f low theory, that of u at the centre 

line, O) / and that of kinetic energy of the disturbance 

gg(X), in Figure 6.16 for X = 0.001 and R = 16000. It is 
observed that for any w, g^ (X, O) is maximum and is positive 
for the widest range of frequencies. The growth rates based 
on ijj and u at the centre line of the channel, on the energy 
E of the disturbance and cn the parallel flow theory for 
different Reynolds numbers and frequencies, and at different 
X are tabulated in Table B4 . Figures 6.17 through 6.21 show 
variation of these different growth rates with oj at X = 

0.001, 0.002, 0.004, 0.006 and 0.008 respectively for diff- 
erent Reynolds numbers. One can see that the above conclu- 
sicn is true for all ccmbinations of X and R- 

The neutral points correspohding to various distur- 
bance properties can be easily found from Figures 6.17 
through 6,21 as they correspond to growth rate being zero. 

The peaks of the growth rate curves for a disturbance 




# 
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property at any X may be joined; the intersection of such a 
curve with the w =aocis determines the critical frequency. 

The critical Reynolds number for the same disturbance 
property at any X is then obtained by plotting its maximum 
growth rate at that X against R (see Figure 6.22) and 
reading the intercept of the curve with the R=axis. The 
data corresponding to the neutral and critical points 
obtained by this procedure are tabulated in Table b5. 

Figure 6.23 shows the different neutral curves, at different 
X, based an (i) g^ at the centre line of the channel, (ii) 
growth of energy of the disturbance and (iii) the parallel- 
flow theory. Also shown is the neutral curve based on 
at the centre line of the channel only at X =0,001 since 
it is little different from that corresponding to (X, O), 
It is observed that the neutral curves are different for 
different flow quantities. However, the difference between 
the neutral curves corresponding to different disturbance 
properties decreases as X increases. The ncnparallel 
effects make the flow unstable at lower Reynolds number and 
for a wider range of frequencies compared to those obtained 
from the parallel flew theory. The actual amount of such 
an effect depends on the choice of growth rate used for 
determining the neutral curve. The growth rate (X, O) 
gives the minimum critical Reynolds number at all X. 
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Figure 6.24 shews the variaticn o£ critical 

frequency/ w / and critical Reynolds number, R , as obtained 

CXI the basis of (X, 0), g^ and the parallel flew theory 

against X. It is found that at X = O.OOl, the R^ predicted 

by the parallel flow theory is greater than that corres^ 

pcxiding to g^ (X, 0) by 22.8%, and that corresponding to 

g by 4,5%. At X = 0.008, these differences reduce to 

8.3% and 1,6% respectively. This is to be expected since 

the ncxparallel effects must vanish at large X. Figure 6.25 

exhibits a similar behaviour with the physical coordinate 

X for the critical Reynolds number based on various growth 

rates (note that Figure 6.24 is a log-log plot). At X = 

20.0 the parallel flow theory overpredicts the critical 

Reynolds numbers by 24.6% and 6.2% as compared to those 

corresponding to g (X, O) and g_ , respectively. These 

\p t 

differences reduce to 10,4% and 2,5% at X = 75, 

^ • 2 . Rigid Circular Pipe 

Assuming the developing flow in the pipe to be 
parallel, spatial stability results have been obtained for 
both the Hornbeck and sparrow profiles. However, ncxi- 
parallel effects are considered only for the better 
(Hornbeck) profile. While Huang and Chen [45] do give 
results for temporal stability of the Sparrow profile on 
the basis of parallel flow approximaticx, they cannot be 
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used directly for comparison since the frequency values of 
the neutral disturbances are not available in Huang and Chen 

6.2.1. Developing Velocity Field 

The variation of the basic flow velocity conponent 
U and its gradient in the radial direction obtained for the 
Hombeck and Sparrow profiles at X = 0.0014 (X = 0,003) 
and at X = 0.00616 (X = 0.01) is shown in Figure 6.26. 

The differences in the two velocity profiles and their 
gradients are similar to those for the channel flow (cf. 
Figure 6.1, Secticn 6.1.1). Figure 6.27 exhibits the 
development of the Hombeck 's velocity field as the fluid 
moves downstream. The velocity data for the Hombeck and 
Sparrow velocity profiles at different axial locaticns ave 
shown in Tables A5 and A6, respectively. 

6.2.2. Parallel Flow STA&lhtjy— Central Mode 

It is well known that the central mode is the 
least stable mode for the fully developed flow in a pipe. 
Theijefore, the stability of developing flow was also 
investigated with respect to the central mode, it was 
found that for the developing flow also, the central mode 
tends to approach neutral stability as R increases 
indefinitely and as u tends to zero. Also for this mode 
the two parameters to end R can be combined into cne so 
that their product to^ governs, in an approximate manner, 
the values of Rk . 
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For the least stable central mode the variations 

of Rk and Rk . with uR as the independent variable for 
or oi 

different values of X are shown in Figures 6.28 and 6.29. 
Though Rk^j^ and Rk^^ vary somewhat with different ccmbi“ 
nations of and R such that wr is constant, these 
variations for too small to be represented graphi- 

cally on the scale of Figure 6.28 whereas those for R^<q^ 
are shown by a band I in Figure 6.29. It is observed from 
these figures that while Rk^^. may be related to ojR kiy ^ 
power law the relationship between wR and Rk . is approxi- 

04 - 

mately of exponential type at all values of X. Also at a 

given value of t^R, k . decreases while k increases as X 
^ oi or 

decreases. This implies that though the central mode 
remains stable, it is less stable for the developing flow 
than for the fully developed flow. For instability of the 
developing flow, therefore, we consider only the wall mode 
in the follox-zing sub-sections. 

6.2.3. parallel Flow stability — wall Mode 

Figures 6.30 and 6.31 show the neutral curves, 

( 0 ) vs. R) and (k^^ vs, r), at X = 0.0005,0.001, 0.002, 

0.00323 , 0.0035 , 0.004 , 0.005 and 0.00616 for the Hombeck 

* 

profile and at X = 0.002, 0.003, 0.004, ,0.006, 0,007, 0.009 
and 0.01 for the sparrow profile. The data corresponding 
to these neutral curves are also given in Tables Cl and C2 
respectively. Our results in Figure 6.31 for the Sparrov/ 
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profile do match very well with the (k -- R) curves of 

or 

Huang and Chen [45] ; the comparison of (w - R) curves, 
however, is not possible because Huang and Chen carried out 
the temporal stability analysis and did not report any 
frequency values. Comparison of the neutral curves for 

'k 

both the velocity profiles at X = 0.00323 (X = 0.006) and 
at X = 0*00616 (x = 0.01) shows that the instability of the 
developing flow in a pipe is set in at lower Reynolds 
number and the area of unstable region enconpassed by the 
neutral curves is larger for the Hombeck profile, similar 
results were obtained for the stability of developing flow 
in a channel and confirm the fact that the flow instability 
characteristics are very sensitive to the mainflow velocity 
field. It is also observed that for both the velocity 
profiles the slope of the lower branch of the neutral 
curves is quite large at small x and that it decreases as 
X increases. While the area of unstable region, critical 
frequency and wavenumber decrease with increasing X, the 
critical Reynolds number first decreases upto a certain X 
and then increases as X continues to increase. 

Figure 6.32 shows the variation of critical 
Reynolds number, critical wavenumber and critical frequency 
with, X for the Hombeck and sparrow profiles. It also 
contains seme experimental results for obtained by 
Sarpkaya [48]. It is observed that the minimum critical 
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Reynolds number is about 11700 for the Hombeck profile, 

while it is about 19800 for Sparrow profile? the vs. X 

curve for the Hombeck profile lying well below that for 

the Sparrow profile (Huang and Chen [45] report a value of 

19900 for the minimum It is, however, felt that the 

results reported here are more accurate (cf. section 5.4, 

since Ay = 0.0025 in the present analysis while Huang 

and Chen took AY = 0.005 upto X = 0.0014 and increased it 

to 0.01 at X = 0.00616). As X increases, the critical 

Reynolds number passes through a minima for both the 

velocity profiles and this minima occurs in the vicinity 

of X = 0.0035. However, the portion of the curve near the 

minimum R is flatter for the Hombeck profile. The ratio 
c 

of the critical Reynolds number for the Sparrow and 
Hombeck profiles increases frcm approximately 1,69 at 
X = 0.0035 to about 1.84 at X = O.OOl and to about 2.21 at 
X = 0.00616, This implies that the Sparrow profile tends 
to becone stable quite rapidly in ccxnparison to the Hombeck 
profile. In other words, the developing flow in a pipe 
with the Hombeck profile will exhibit instability in a 
greater length of the pipe. This agrees partially with the 
experimental observations of sarpkaya [48] that the flow 
remains unstable further downstream than that predicted 
by Huang and Chen [45], We also note frcm Figure 6.32 
that sarpkaya* s experimental results for R^, though well 
below the R^ vs. X curve for the Hombeck profile, are 
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still closer to it than to the R vs. X curve for the 

c 

Sparrow profile. Since the flow in the developing region 
is actually nonparallel, it appears that the parallel flow 
assumption is not quite valid at least in the near entry 
region, sarpkaya found the minimum critical Reynolds number 
to be about 3800 for X in the range 0.012 X 0.02. 
However, as sarpkaya himself noted, his critical Reynolds 
number may be low due to the superpositicn of seme nen- 
axisymmetric disturbances on axisymmetric disturbances as 
well as due to a higher initial disturbance level than that 
warranted by the linear theory. In terms of physical 
length x in the entry region, the first instability of the 
flow is theoretically found to occur at X - 37 for the 
Hombeck profile and at X - 64 for the sparrow profile, 
while sarpkaya' s experimental results yield a range of 
45 X 4 75 for the sane. This is evident from Figure 6.33,. 
which further shows that the Hombeck profile gives results 
closer .to the experimental ones. 

One also notes from Figure 6.32 that the critical 
frequency and critical wavenumber for the Hombeck profile 
are larger than the corresponding values for the Sparrow 
profile at all x and that with increasing X both decrease 
first sharply and then gradually. The dimensionless phase 
velocity, w critical point increases from 

0,4 2 and 0,37 at X = 0.001 for the Hombeck and Sparrow 
profiles, respectively, to 0.49 and 0.45 at X = 0.00616. 
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However, these values are much lower than that for the 
fully developed flow for which Garg and Rouleau [40 ] found 
the phase velocity to be slightly less than 2.0. One may, 
therefore, surmise that as the flow develops downstream 
the critical frequency and wavenumber change so that their 
ratio increases towards the value for the fully developed 
flow. Unfortunately, Sarpkaya [48] has not provided 
experimental results for critical frequency corresponding 
to an axisymmetric disturbance. However, he does tabulate 
the experimental values of for non=axisymmetric distur- 
bances and the theoretical values of obtained by Huang 

[102] for the Sparrow profile, it appears from the 
comparative trend of theoretical values of for the 
axisymmetric and non -axisymmetric disturbances for the 
Sparrow profile that the critical frequencies found here 
for the Hombeck profile should compare very well with the 
experimental values for the axisymmetric disturbance. 

sarpkaya [48] points out that his experimentally 
observed velocity profile agrees with the Sparrow profile 
to within 5% whereas Figure 6.26 shows that the Hombeck 
profile also agrees with the Sparrow profile to within 5% 
at sections closer to the entry section. It would have 
been interesting to compare the experimental and Hombeck 
velocity profiles, specially when it is found that the 
Hombeck profile gives stability characteristics that are 
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closer to the experimental ones. Unf ortunately^ however, 
Sarpkaya's velocity profiles, are not available. 

Corresponding to the neutral point for w = l.o 
and R = 11779 at X = 0.0035, Figure 6.34 displays the real 
and imaginary parts of the eigenfunctions for the x and Y 
components of the disturbance velocity. These have been 
normalised with respect to their respective maximum 
absolute values; the ratio of vS to !u| _ being 0.2432. 
It is obvious that the boundary conditions (3.17) are 
satisfied. 

Figure 6.35 shows a representative variation of 
vs. u at X = 0.001 for R = 16000, at X = 0.002 for R = 
14000 and 16000 in addition to such variaticns for R = 

18000 at X = 0,001, 0.002, 0.00323, 0.004 and 0.005. It 
is observed that the k^^ vs. to plot is almost a straight 
line and the slope of the straight lines varies little 
with X. It is also observed that the slope of these lines 
at any X increases slowly with Reynolds number. This 
behaviour helps a great deal in determining the starting 
value of k close enough to the actual one so that the 
convergence is achieved in a few iterations. 

6.2.4. Effects of Nonparallelism of the Flow 

As noted above the oily experimental work on the 
stability of developing flow in a rigid circular pipe is 
due to Sarpkaya [48]. Though he mentions that the 
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streamwise component of the disturbance velocity was 

measured at different radii, he neither reported the magni= 

tude nor the radii at which measurements were made . Since 

these details could not be obtained [95] , it was decided t; 

compute all the growth rates mentioned earlier. Neutral 

curves were found on the basis of q_ and the values of g 

"E 

and g^ at the pipe axis since g^ and g^ are functions of 
Y also. H'Cwever, it is clear from eqns. (4.4) and (4,6), 
using appropriate boundary conditions at Y = O, that 


g (X, O) = g (X, O) . 
u \p 

Growth rates based on u, and E were obtained at 
X = 0.0005,0.001, 0.002, 0.0035, 0.005 and 0.007. Figure 
6.36 shows the growth rate based on as a function of Y 
for different combination of X, R and u. It is observed 
that the maximum growth rate occurs at the pipe axis and 
the variation of g with Y is similar to that in the case 
of developing flow in a channel (cf. Figure 6.15). Thus 
remarks made in Section 6.1.4 for Figure 6,15 hold here as 


well. 

Figures 6.37 through 6.42 show the variation of 
g (X, O) and g„(X, O) with w at X = 0.0005, 0.001, 0.002, 

}p t. 

0.0035, 0.005 and 0,007, respectively, for different 
Reynolds numbers. Just as for channel flow, it is obsearved 
that for any w at a given X and R/ 9^ (X, 0) is greater 

than g_ and is positive for the widest range of frequencies. 

' 'E' 

One can also obtain the data corresponding to the neutral 
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and critical points in this case following the method 
outlined in section 6.1.4 for the channel flow. Figure 6.4 3 
helps in finding the critical Reynolds number in the 
present case. The neutral and critical points obtained in 
this way are tabulated in Table c4 . Also Figure 6.44 shows 
the different neutral curves at various based on (i) 
g (X/ O), (ii) g^(X) and (iii) the parallel flow theory. 

Tp ^ 

Similar to the channel flow case, it is observed • that in 
comparison to the results for the parallel flow theory, the 
nonparallel effects make the flow unstable at lower 
Reynolds number and for a wider range of frequencies. The 
actual amount of such an effect depends on the choice of 
the growth rate used for determining the neutral curve. 

The growth rate g (X, 0) gives the minimum critical 
Reynolds number at all X- 

Figure 6.45 shows the variation of the critical 
frequency, and the critical Reynolds number, R^, as 

obtained on the basis of g, (X, 0), g (X), and the parallel 
flow theory against X. It is observed that the variation 
of with X is similar f or various disturbance properties; 

corresponding to g^ (X, 0) being maximum and that 
corresponding to the parallel flow theory being minimum. 
Also, all the R^ vs. X curves pass through a minima. The 
minimum critical Reynolds number corresponding to g (X, 0), 
to g (X), and to the parallel flow theory are 9700 at 

ill 

X = 0.00325, 11000 at X = 0.0035, and 11700 at X = 0.0035, 
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respectively, in comparison to the results based on 

(X/ O) and g„(X)/ the parallel flow theory overpredicts 
the critical Reynolds number by 29,&/o and 3.7% respectively 
at X = 0.0005, by 20. C% and 6.4% respectively at X = 0.0035 
and by 26.5% and 12.0% respectively at X = 0.007, This 
implies that the vs. X curves obtained on the basis of 
nonparallel flow theory are flatter than those corresponding 
to the parallel theory; R^ does increase beyend X = 

0.0035 but not so sharply as for the parallel flow theory. 
Physically, it means that the actual developing flow is 
unstable over a larger inlet length of the pipe than its 
parallel flow approximate. The first instability of the 
flow on the basis of g, (X, O), g (X), and the parallel flow 
theory is found to occur af X = 33, 38, and 37 respectively 
(see Figure 6.46). The R^ vs. X or X curve obtained on the 
basis of g (x, 0) is closest to the experimental data of 
sarpkaya [48] . 
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Chapter 7 
' CONCLUSIONS 

Spatia.1 stability characteristics of the developing 
flow in a two dimensional channel and in a rigid circular 
pipe have been obtained for the actual nonparallel flow as 
well as for its parallel flow approximate. The developing 
flow velocity field has been determined by the finite 
difference method, considering the developing flow in the 
channel as a parallel flow# it is found that the symmetric 
disturbances are more unstable than antisymmetric distur- 
bances. Thus, for most of the analysis, the disturbance is 
assumed to be symmetric for channel flow, and axisymmetric 
for pipe flow. Further, only the wall modes have been 
considered in case of pipe flow stability since they exhibit 
instability in the developing flow regicn. The central mode, 
on the contrary, remains stable. This mode is, however, 
fomd to be less stable for the developing flow than for 

the fully developed flow. 

For the parallel flow analysis, results have 
also been compared with those for the velocity profile 
given by Sparrow et al.'s linearization method. This 
comparison reveals that even though the two velocity 
profiles differ by less than 5 % throughout the field of 
flow in the near-entry region, their stability characteristics 
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are quite different. Also, the comparative 
stability characteristics obtained for both 


^ of the 
^tries of 


flow on the basis of parallel and nonparallel 
that the results are significantly different. 


th 


Tories shows 


non 

parallel effects introduce the following main 
(i) The growth rate beccmes a function of the 


as well as the transverse coordinate. The 


®ren ces: 
^treamwise 


9towth 


rate of the disturbance stream function^ , , 

^ich is 

maximum at Y = 0, diminishes gradually 

Uniformly 

towards the boundary layer edge and change 

® rapidly 

in the boundary layer region; the amount 

variation 

of the growth rate in this region is a f pj-, _ 

^tion of 

the axial location and frequency of the 

®turbance. 

(ii) The growth rate also becones a function qj 

disturbance property involved, that is, tha ^ 

^ stream 

function, velocity components, and kinetic 
density of the disturbance have different 
rates. Thus at any point in the flow reg^^ 
disturbance property may be growing while 
be decaying. This implies that there are 
neutral curves and hence different criticej_ 
number and frequency different disturh^^^ 

properties. 

(iii) The actual developing £1°^ becones unstshi, ^ 

wider range of frequencies as compared to 

parallel flow approximate; the actual of 


such an effect depends cn the growth rate used for 
determining the neutral curve, it is found that 
g^ (X, 0), the growth rate of the disturbance stream 
function at Y = 0, (a) is the maximum of all growth 
rates at any frequency for all r and X, (b) remains 
positive for the widest range of frequencies, and 
(c) gives the minimum critical Reynolds number at 
all X. 

(iv) Like the growth rate, the disturbance wavenumber also 
becomes a function of the streamwise and transverse 
coordinates as well as a function of the disturbance 
property. However, the wavenumber obtained from the 
parallel flow theory is little affected by the non-< 
parallel effects. 

It is found that the step-by~step integration of 
the continuity and mcmentum equations separately instead 
of a single differential equation for the disturbance is 
more economical as well as accurate. Also, selective 
application of the Gram~Schmidt orthonormalization procedure 
is found to be the best way to control the parasitic 
errors during numerical integration. 

The following conclusions apply to the two 
geonetries of flow separately. 
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7.1. parallel Plate channel 

When compared with the parallel flow stability 
results for symmetric disturbances, antisymmetric disturb 
bances are found to be more stable at all x except for some 
eigenstates that lie close to the lower branch of the 
neutral curve and correspond to Reynolds numbers much greater 
than the critical value; the critical Reynolds number being 
lower at all X for the symmetric disturbances. The critical 
Reynolds number for the symmetric disturbance decreases 
continuously with X while for the antisymmetric disturbance 
it first decreases to about 12100 at X = 0.005 and then 
increases rapidly so that the fully developed flow seems to 
be stable to all antisymmetric disturbances. 

For the symmetric disturbances, the developing 
flow is found to be more stable than the fully developed 
flow. The critical Reynolds number, wavenumber and 
frequency decrease with increasing X and approach asympto- 
tically the corresponding values of 3848.1, 1.0198 and 
0.4O369 respectively for the fully developed flow. The 
comparative study of the stability behaviour of the B-O 
and the Sparrow velocity profiles for the developing flow 
considering it as a parallel flow shows that the critical 
Reynolds numbers are same for both the velocity profiles 
for X $ 0.084 (or X 05= 440) but for X < 0.084, the critical 
Reynolds number for the B-0 profile is less than that for 



the Sparrow profile owing to the difference in the two 

velocity profiles. The critical Reynolds number at X = 60 

for the B=0 profile is about 10900, which is almost half of 

that for the Sparrow profile. This large difference in 

the critical Reynolds numbers implies that the shape of 

velocity profile is an important factor for deciding the 

stability behaviour of the flow. In the absence of any 

experimental data on the stability characteristics of the 

developing flow in a channel, it is difficult to determine 

0 

which velocity profile predicts the stability of the real 
flow more accurately. However, since the superiority of 
the B-0 profile over the Sparrow profile has been accepted 
by several research workers [88=93] , one has the intuitive 
feeling that the stability results for the B-0 profile 
should be closer to the actual cnes. 

The nonparallel theory predicts critical Reynolds 
numbers which are lower than those obtained from the 
parallel flow theory. The latter predicts a value which is 
higher than that corresponding to 9^ (X/ O) and to for 
the nonparallel flow by 22,8% and 4,5% respectively at 
X = 0.001, These differences reduce to 8.3% and 1,6%, 
respectively at X => 0,008. Moreover, the neutral curves 
based on different growth rates tend to merge into the 
neutral curve for the parallel flow theory as X increases, 
implying thereby that the nonparallel effects vanish for 
large X as they should. 
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• 2 ♦ Circular Pipe 

Like the fully developed flow, the developing flow 
is stable to the central mode. However, unlike the fully 
developed flow, it is unstable to the wall mode. At all 
axial locations the critical frequency and wavenumber for 
the wall mode obtained cn the parallel flow approximation 
for the Hombeck profile have been found to be greater than 
those for the Sparrow profile. However, reverse is true for 
the critical Reynolds number. The critical Reynolds numbers 
at X = 0.0035 for the Hornbeck and Sparrow velocity profiles 
are respectively 11700 and 19800. As one moves downstream 
or upstream of this point, the vs. X curve for the 
Sparrow profile rises more rapidly in conparison to that 
for the Hornbeck profile implying thereby that the Hornbeck 
profile is unstable over a larger inlet length of the pipe. 

The actual nonparallel flow in the developing flow 

region of the pipe yields a lower critical Reynolds number. 

The minimum critical Reynolds number corresponding to 

g, (X, 0) and to g„(X) for the Hombeck profile are 9700 at 

yj E 

X = 0.00325 and 11000 at X = 0.0035, respectively. The 
parallel flow theory, in comparison to the results based 
cn g^ (X, O) and gg(X), overpredicts the critical Reynolds 
number by 29.8% and 3.7% respectively at X = 0.0005, by 
20% and 6.4% at X = 0.0035, and by 26.5% and 12.0% 
respectively at X = 0.007. The R^ vs. X curves obtained 
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on the basis of nonparallel flow theory are flatter than 
those obtained from the parallel flow theory. Thus the 
actual developing flow remains unstable over a larger inlet 
length of the pipe than its parallel flow approximate. 

The first instability of the flow, on the basis of 
g^ (X, 0), and the parallel flow theory, is found to 

occur at X =; 33, 38 and 37 respectively. The vs, X or 
X curve obtained on the basis of g (X, O) is closest to 
the experimental data of sarpkaya [48] . 
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APPENDIX A 

DEVELOPING VELOCITY FIELD IN CHANNEL AND PIPE 


The Finite Difference Method 

The finite difference method for determining the 
velocity field in the flow through a channel was developed 
by Bodoia and Osterle [ 20 ] . Later Hornbeck I 47 ] extended it 
to the pipe flow case. They took a boundary layer model for 
the flow in the entrance region and replaced the derivatives 
by finite difference relations. This results in a set of 
simultaneous algebraic equations. The solution of these 
equations determines the velocity field in the flow region. 
The details of the method for the channel and pipe flow are 
given below. 

For an incompressible , Newtonian fluid, the 
boundary layer equations for the basic flow, using the ncn- 
dimensional variables defined in eqns. (2,9), are 


and 
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where m = 0 for the channel flow and m = 1 for the pipe 
flow. The new nondimensional variables used here are 


defined as 
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and 


X = X/R , 

P - P 

« _ o 



V = VR , 


(A. 3) 


where is the pressure at the inlet to the duct. The 
boundary conditions on U(X/ Y), etc. , a-f«- 


and 
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at Y = 1 for all X b > 

at Y = 0 for all X o (A. 4) 

at X = 0 for all Y . 


Equations (A.l) and (A. 2) with (A. 4) are solved 
using the finite difference method. A grid as shown in 
Figure A-1 is developed. At any grid point (j, k), using 
implicit finite difference representation, one can write the 
finite difference form of the eqn. (A. 2)^ as 
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(A. 5) 

for k = l(l)n, where (n+2) is the number of mesh points 
along Y at any cross-section j. Also, the finite difference 
form of eqn* (A.l) is. 
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for k = l(l)n. Equations similar to (A. 5) and (A, 6) can be 
written for k = 0 in case of channel flow invoking the 
symmetry conditions at the centre line and using the fact 
that m = 0 in this case. However, in the case of pipe flow, 
term involving Y in denominator of eqns. (A.l) and (A. 2) 
need to be evaluated at Y =! O by L'Hospital rule, when eqn. 
(A. 6) is written for k = 0 -o n for channel flow or 1 to n 
for the pipe flow and the resulting equations, including the 
one obtained from L'Hospital rule in case of pipe flow cnly, 
are summed up, we get the following simple equation 
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The finite difference form of the resulting equation frcm 
the application of L'Hospital rule to eqn. (A. 2) is 
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Equations (A. 5) with k = 0 also and (A. 7) for the channel 
flow or (A.5), (A. 7) and (A. 8) for the pipe flow give (n+2) 
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equations in (n4-2) unknowns viz., U's at k = 0(1 )n and p at 
the cross-section j+1. These equations form a matrix which 
is of tridiagonal type except for the row corresponding to 
eqn. (A. 7). For solution, this matrix may be partitioned to 
take advantage of the tridiagonality of the matrix. The Y 
component of velocity can be obtained from eqn. (A. 6). 

It is observed that this method is a marching 
method and therefore requires the values of U, v and P at all 
the grid points at the entry section. Since the discreti- 
zation of the duct cross-section, taking , = 1 for k = 
0(l)n, leads to a reduction in the volumetric flow rate 
because of the no slip condition at the wall, a form factor 
such as 




- AY/2) 

for channel flow , 

(A. 9a) 

or 

Uo,k = 

- 1) 

for pipe flow , 

(A. 9b) 


For k = 0(1 )n was used [ll3]. 

The velocity profiles for the developing flow in 
the channel and in the pipe are tabulated in Tables A3 and 
A5 and shown in Figures 6.2 and 6.27, respectively. 

A. 2. Linearization Method 

The linearization method was first proposed by 
Langhaar [77] and was used later by Sparrow et al. [l9] 
to calculate the velocity distributicaa in the developing 
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flow through a channel and a pipe. They recasted the momentum 
eqn , (A. 2) into the form 


e*(X)u3^ 


9_U 

9x 


A(X) + ( 




m 

Y 9Y^ 
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(A. 10) 


where e (x) is a weighting function of X to be determined and 
A(X) is another undetermined function which includes the 
pressure gradient as well as the residual of inertia terms. 
Finding A(x) by integrating eqn. (A. 10) over the cross- 
secticn while using the continuity equation and introducing 
the stretched streamwise coordinate x , defined by 

— k k 

dx = e d X , (A. 11) 


we get the following monentum equation 
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The boundary conditions are 

U = 0 at Y = l, 3-2 = 0 at Y = 0, 

o Y 

(A.13) 

k 

and u = 1 at X = 0 . 

By expressing the velocity U at any point in the 
inlet region as the sum of fully developed flow velocity 
and a difference velocity. Sparrow et al. get a differential 
equation which has an infinite convergent series solution; 
the constant of the series being determined by using the 
orthogonality conditions. Thus they arrive at the 
following solution. 
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where, A^ are the roots of 

J^( A^) » 0.5 )v. J^(A^) = 0. (.Y.17) 

Vf 

The weighting function e for both the flows is given by 

/ (2U - 1.5 U^)(aU/9X*) y”' dY 

* o 

e = j — . (A. 18) 

[9U/)y]y^^ + / (9U/3Y)^ Y^ dY 

o 

The first fifty values of A for the channel and pipe flow are 
given in Table Al while values of X corresponding to a few 

'k 

values of X for both the flows are given in Table A2, The 
velocity profiles for the developing flow in the channel 
and in the pipe are tabulated in Tables A4 and A6 and 
shown in Figures 6.2 and 6.21 , respectively. 
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faole A?. ReiattonshiD between X and X 
for channel and pipe 
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ftpnen.iix B 

'Afif.fcS OF STABII.ITY RKSUfjTS FOR THE CHAMRFt. FLOW 


’risvie .ji The neutral stability results for the B-0 
profile and antisv.nfnetric dlsturbanres 
Cnaraliei flow tneoryi 
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2 , n U 2 4 3 

1 1 - s /f 

» ? i/ 1 

/I 1 C 3 i a 

t .oyouu 

72265’" 

22000, 

1.00487 

2.5b269 

# ' H ^ 

/;?u‘!o , 

1 , 1 'Jb7R 

2 . B 1 7 d 3 

2 6 u 00 , 

f|.7b099 

1 ,99116 

• *>) u ^ 

u ^*1/ » 

1 a ^ 19 

3.12167 

30000, 

0,64362 

1.71924 

m 0 U 1 

i*'ooo. 

U t r.i A 7 ^ 

3 .IbHbb 

34000. 

().5b245 

1 .52462 

. h 3 

a 1 ^ ! bH5R 

.3 ,1 64 65 




.OUS 

4.12 0,,). 

1 .7 49 40 

23200, 

1 ,02868 

2. 66473 

> 1! Uh 

*r,U Oy , 

6u2 u7 

1 ,S29bl) 

26000, 

1.02937 

2.65701 

. h ye 

1 oyeo . 

0, 7dl 4t 

1 ,88797 

21 400, 

0,66771 

1 . 5930 " 

.hub 

21'* 00, 

9 * H 8 1 tr/ 

2 .1o755 

25000. 

0.53090 

1.29969 


2 So 0(1 , 

0, ”91 V7 

2.25244 

19900. 

0,81639 

1.91604 

.hub 

iOgOU. 

fi,0-2S()0 

1 .06202 

30000. 

0.86383 

2.2^927 

. 0 U7 

20221. 

0,700UO 

1 .66934 

20400. 

0.65427 

1154542 

• 0U7 

2040U, 

0,7 23 51 

1 ,72192 

240bi), 

0,49537 

1.1 7275 

,0U7 

2'H»0U, 

/7 

i .89012 

2800U, 

0,40972 

0 , 97.107 

,0 07 

2Si./no. 

0 , 7 o? bo 

1 ,99921 
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.Ouo 

2Sy3u, 

0 , .i<jO y'o 

1 .1 0422 

25b00. 

0,43513 

0,99207 

,0 0'^ 


o.stmtii 

1 ,164 9 0 

30000. 

0.34101 

0,77270 

.0 09 

.30u()u , 

0.5iSo3 

1 .23713 

34000, 

0,29263 

t ^ t> 7 C) 7 

.0 00 

i4u0y , 

0,4 07 IS 

1 .21941 
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.010 


n, 3 6000 

0,82778 

34000, 

0.32179 

0,728S8 

.010 

i^OOlla 

0 , 3 46 

0, 84Bu4 

38000. 

0.27038 

0,60378 

. 0 1 0 

3H0nu, 

0 » ..i on 34 

0.86636 

42000, 

0.23314 

0,53332 

. 0 i 0 

42U0U. 

0,. 35434 

0.84997 
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Idhii? C ,1 ■'* d vonuffiDP r ».ind rates tor the HorntecK 

r-roti,le ana axi symmetric disturbances 
iparailel and nonparallel flow) 


xs 0,0005 


21 i)U*>,0 
710^0,0 
71000.0 
2 aU0U.() 

71000.0 

71 000.0 
7lOou,i! 

71000.0 
7;iOu!?.o 
?.lOuO,n 

70000.0 
70000 , 0 
7UUofi.O 
7 0 0 0 0 . 0 
70000,0 
70000,0 
700oU,0 

2 U 0 0 0 , 0 

70000.0 

20000 . 0 
20000.0 
I y 100.0 
lyduo.o 
19300,0 
,19300.0 
19300,0 

19300.0 

laoou.o 

18000.0 
18000,0 
1 8000,0 
18000,0 
18000,0 
Id 000,0 
18000,0 
IbOOO, 0 
10000,0 
1 b 0 0 0 , 0 
1 0 0 U 0 , 0 

1 b 0 0 0,0 
) b 00 0,0 
IbOwU. 0 
1 b 0 y 0 , 0 


Gj 


2.00 
2,20 

2.40 
2,bO 
2.80 

3.00 

1.10 
.3,15 
3.30 

3.40 

7 . 1 0 
2,70 

2.40 

2. b0 

2.7 0 
2.80 

2.90 

3.00 
3,20 

3.40 

3.50 
2.7 0 
2,80 

2.90 

3. 00 

3.10 

3.50 

2.40 

2.50 

2.50 
2.80 

3.00 
3.30 

3.50 

3.50 

2.50 
2,80 
3.00 
3,20 

3.40 

3, U0 
3,20 
3.40 


nr 

57*141017 
5,581892 
5,019923 
6,455596 
5,889414 
7.321lo9 
7,535292 

7,543555 

7,954985 
,1 /8480 
,3 31597 
5,549901 
5.984455 

5,415548 
,532051 
,845946 
7,051373 
7,2/5321 
7.701811 
125290 
337529 
,502/10 
5,815322 
7,029521 
7.242248 
7,454530 
8.298844 
5.909911 
5,927881 
5,334539 

5.75/473 

7.178487 
7.805959 
8,223/23 
y. 4313/7 
5,245484 
5,551085 
7.0/4019 
7,485445 
7.895345 
7,018517 
7.425570 
7,831113 


8 , 

5, 


5, 

6, 


8 , 
8. 
5 , 


oi 

'07057737’ 

0,024511 
0.000513 
■'0,013521 
•0.01/052 
•0,009556 
-0,001514 
U. 003611 
0,0234/6 
0.040552 
0,060416 
0,044383 
0,018916 
0,002841 
•0,001447 
•0,003196 
•0. 002310 
0.001243 
0.015651 
0,043478 
0,061283 
0,010379 
0,007729 
0,007603 
0,010057 
0.015188 
0,063143 
0,061525 
0,041791 
0,041794 
0.030972 
0.029557 
0.045024 
0,070004 
0,085029 
0.090524 
0,075425 
0,068565 
a;070594 
0,081602 
0.093122 
0,092028 
0.099379 


0 

E 

•0 7059785 “ 
-0,024152 
0,002579 
0.019966 
0.027238 
0.023717 
0.017729 
0,013780 
-0.002978 
-0.017857 
-0.061303 
-0.044410 
-0.016189 
0.003054 
0.009069 
'0.012674 
0.013787 
0.012174 
0,001051 
-0.021538 
-0,036975 
-0,003102 
0.001242 
0.003327 
0,002842 
-0.000118 
-0,039489 
-0.059702 
-0,025352 
-0.035996 
-0,022755 
-0,018225 
-0.027868 
-0,047394 
-0,061160 
-0,086653 
-0,070390 
-0,058088 
-0.055312 
-0,063093 
-0.083472 
-0,078632 
-0.082070 


q CX,0) 
y' 

■o7oTl574 

•0,004074 
0,050815 
0.06/9/2 
0,074976 
0,071246 
0.065205 
0.061277 
0,0445/3 
0,029817 
•0,009984 
0.005976 
0.035133 
0.054159 
0,060028 
0,063483 
0,064448 
0.062716 
0.051428 
0.028977 
0.013722 
0.050207 
0,054447 
0.055356 
0,055727 
0.052634 
0.0U261 
-0,001055 
-0.001396 
0.021490 
0.035427 
0,039701 
0.029543 
0.009925 
-0.003775 
-0. 018676 
-0.002544 
0,0092/9 
0.U10618 
0,003465 
-0,010230 
-0,005845 
-U, 0097 27 
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12f'U0. U 

2,00 

4.604580 

0,0652/9 

-0,053470 

-0,01 3668 

lii'un.y 

2,20 

5,0100/3 

0,057380 

-0,044302 

-0,004357 

.woyf*,') 

2,40 

5,414743 

0,060694 

-0,045900 

-0,005908 

12SUU,(» 

2,00 

4,62/859 

0.053550 

-0,042042 

-0,0041/3 

12bt>0,U 

2.20 

5, 036971 

0,047260 

-0 , 034406 

0. 003580 

i2bui:,o 

2.40 

5.443427 

0,05267.3 

-0,038043 

-0,000012 

1, 4(jU'UU 

1.70 

4,062835 

0.048154 

-0.038769 

-0, 0061 25 

1 4nt*!».u 

2,0 0 

4,694209 

0.023655 

-0,012833 

0,020145 

3 ‘tOUU.O 

2,20 

5.110953 

0.022516 

-0,0101 62 

0.022892 


2. iO 

5,318259 

0, 026835 

-0.013565 

0,019519 

1 "J! (} V ' f * 0 

2.bO 

5. 730811 

0,045755 

-0,030344 

0. 002838 

1 'iOU(> , 1.) 

2,bb 

5,833515 

0,052698 

-0.036689 

-0.003465 

laouo.o 

1.80 

4,311146 

0,01 8236 

-0,008941 

0.021240 


1 ,yo 

4,523937 

0,011279 

-0.001438 

0.U28B24 

.1,3«>UU,y 

2.00 

4, / 3584.3 

0,00/515 

0.002984 

0,033302 

1‘30m(>.0 

2, if) 

4,946938 

0,007065 

0,004200 

0,034557 

ibouo.o 

2,20 

5,157271 

0,010044 

0.002085 

0,032478 

IbOuO.U 

2.iO 

5.366865 

0,016561 

-0,003468 

0,026966 

1 oouu , 0 

1.50 

3.694/34 

0,037709 

-0,029828 

-0,002338 

IbOOO.O 

1.7 0 

4,130414 

0,010668 

-0,002203 

0.025619 

loOOCJ.O 

1,80 

4.346450 

0.001763 

0.007167 

U,U35106 

IbOuO.O 

1 .yo 

4.561490 

-0,003858 

0,013408 

0,041392 

i&uyo,!,) 

2.00 

4.775632 

-0.006059 

0,016315 

0.044346 

InOvJO.y 

2,10 

4.988940 

-0,004714 

0,015776 

0.043848 

Ibouu.o 

2,20 

5,201449 

0,000298 

0,011669 

0,039785 

ibouu , 0 

2,40 

5.624088 

0,021790 

-0,007715 

0.020553 

160 Ob, a 

2.55 

5.938893 

0.048402 

-0,032472 

-0,003985 

18000,0 

l.iO 

3,298445 

0,041036 

-0.034224 

-0,010735 

18000,0 

1,40 

3.524485 

0,021453 

-0,014591 

0.009150 

1 8000.0 

1.50 

3. 748669 

0,004829 

0,002240 

0,026168 

18000,0 

1,60 

3,9/1315 

-0,008601 

0.016015 

0,040077 

I aouu.o 

1,80 

4,412869 

-0.025103 

0.033564 

0,05/788 

18000,0 

2,00 

4,850307 

-0.026714 

0,036638 

0,060963 

1 8000,0 

2,20 

5,284087 

-0,012285 

0.024062 

0.048521 

1 8000,0 

2,30 

5,499611 

0,001282 

0.0U570 

0.0361 41 

1 8000,0 

2,40 

5.714176 

0’019267 

-0,005234 

0,019498 

1 8000,0 

2,50 

5,92/717 

0.041809 

-0,026479 

-0,001523 
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X,= 0,0020 


K 

U) 


K 

a 

q {X,0) 



or 

01 

t 

Y 

ToouoCi*'’'™''”' 


~i7l022iB““ 

0.030825 

-0,6^018 

-0.010727 

1, 0 il y ' 1 ® il 

1 , bO 

i t 6 1 U 6 0 2 

0.034526 

-0.0248/5 

-0,0042 IB 

1 0 0 u 0 , 0 

.1 . 70 

3. 712H34 

0,03/602 

-0,027341 

-0.006624 

1 U*>Ui! , (,) 

1 # ‘^■S 

3.225a29 

0,028401 

-0,010821 

-0,0005/0 

lu^o^,y 

1,55 

3. 431503 

0.026849 

-0,017721 

0,001669 

1 0 b 0 , i„i 

i # 65 

3.635082 

0.020107 

^Q^Q1948R 

0.000020 

12000,0 

1.20 

2,75iy73 

0.022036 

-0,01 4881 

0.001211 

i20uo.n 

i,iO 

2.36(5502 

0,013319 

-0,005926 

0,010343 

1 0 

1,40 

.3, 1 /yoy4 

0.008545 

-0,000789 

0,015624 

1 20 oo,u 

l,b(i 

3, 340004 

0, 008007 

0,000235 

0,016769 

1 iUlHI.fJ 

1,00 

3,5yy4i2 

0,011055 

-0,003110 

0.01 3531 

12000,0 

1 . bO 

4,014116 

0.034249 

-0.023821 

-0.006971 

1 

l.iO 

2,yyby5o 

0.000273 

0,006716 

0.021548 

1 iOUO.O 

1 .4:0 

3,214232 

-0,002727 

0,010131 

0.025089 

1 Ji 

1 .50 

3,42/806 

-0.001103 

0,000039 

0,024107 

IJOuiKV 

1 .bO 

3.530824 

0,005415 

0.003174 

0,018343 

1 40in.l*0 

1 , u 0 

2.363305 

0,010548 

-0,013530 

-0,000404 


1 . 10 

2.5H7y59 

0,005450 

0,000655 

0.0149/9 

1^000,0 

1,30 

3.020719 

-0,010556 

0.017224 

0,030842 

14000,0 

1.50 

3,46i5bl 

-0,00/789 

0,015490 

0,020324 

14000,0 

1 » b 0 

3,67/981 

0,001402 

0,006994 

0,020929 

140 00,0 

1. 70 

3.890830 

0,016173 

-0 ,006951 

0,00/096 

14000,0 

1 ,ao 

4,102008 

0,036804 

-0,026611 

-0.012429 

loU00,0 

o.ao 

1 ,038135 

0,031876 

-0,026542 

-0,015668 

loOUO.O 

o.yo 

2.174180 

0,012143 

-0,006901 

0,004201 

loo 00,0 

1,10 

2,635854 

-0,016530 

0,022013 

0,0334o0 

10000,0 

1,30 

3.086000 

-0,026751 

0.032952 

0.044637 

lOOUO.O 

1,40 

3,300251 

-0,023864 

0,030595 

0,042379 

160UO.O 

1.50 

3.520 721 

-0.015176 

0.022555 

0.034439 

10000,0 

1,60 

3.748372 

-0,000355 

0,008512 

0,020509 

10000,0 

1.70 

3.065189 

0,020935 

-0.011852 

0.000284 

1 «000,0 

y . 7 0 

1,725128 

0^033024 

-0,028225 

-0.018901 

IdOOO.O 

U.BO 

1,060517 

0,010935 

-0,006280 

0.003264 

loouu.o 

o.yo 

2.200607 

-0.008102 

0,012765 

0,022496 

liJOUO,U 

1.10 

2,6/0/80 

-0,033251 

0.038304 

0.048319 

laouo, 0 

1.30 

3,130136 

-0,037110 

0,043001 

0,053222 

lauuo.o 

1 . 50 

3,580753 

-0,016132 

0,023331 

0,033759 

1 aouo , 0 

1 . 60 

.3,811043 

0,004653 

0,003408 

0,013972 

lauuo.o 

1.70 

4,031042 

0,032843 

-0.023748 

-0,013005 
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x,= 0, 

0035 



K 


Laj 


k 

Of 

g CX.O) 




Or 

“Oi 

E 

y” 


• y 

1 * U 0 

2, 104855 

0,029007 

-0,022635 

-0, 01 iS'ii 

‘'1 y y 

# 

i m 1 y 

2. 113298 

0,024953 

-0,018290 

-0,00/118 

'4 (1 ij y 

# y 

i <1 4!i y 

2,519549 

0,025344 

-0.018276 

-0.006976 

i u 0 y 0 

* 

1,00 

2. 137360 

0,016261 

-0,010376 

-0,000600 

1 

« i) 

1.20 

2,556903 

0.017185 

-0.010501 

-0,000490 

1. 1501! 

• 

y 0 B 0 

1. 7 4150 3 

0,015081 

-0, 010216 

-0. 0021 32 

1 150*’ 

m 

0*9(1 

1.963117 

0.006273 

-0,001211 

0,00/009 

1 1 ^ y f) 

m 0 

i . 00 

2.182119 

0,002121 

0.003240 

0,011580 

1 i 50<* 

0 0 

1.10 

2,396645 

0,003352 

0.U02414 

0,01 0857 

1 1500 

• 0 

1,20 

2, 612820 

0.010565 

-0. 004281 

0,004253 

■Ubuo 

• 0 

1 , 30 

2,824/21 

0,024296 

-0,017362 

-0.008747 

1 ,10 00 

• 

U.bO 

1,313927 

0,027287 

-0.023083 

-0,016292 

3 iOoO 

• 

u, /o 

1,545111 

0,012499 

-0,008252 

-0,001334 

1 i')UU 

0 Cl 

1.00 

2,222851 

-U,00?3?6 

0.012371 

0.019630 

1 5 0 00 

• 

i. 10 

2,443494 

-0,002537 

0.007987 

. 0.015383 

.1 3'lu!! 


1,30 

2 , b / / 2 1 1 

0,028498 

-0,021743 

-0,014238 

1 40 0 0 

0 Cl 

0,60 

1,328624 

0,018538 . 

-0,014644 

-0,008394 

1 -iOoU 

#C) 

0,80 

1,793 / 73 

-0.006639 

0,010798 

0,017280 

140UU 

• Cl 

1.20 

2,691717 

0,0106/5 

-0,004761 

0,002074 

1 400 0 

# C) 

1.30 

2,909443 

0,033835 

-0,027139 

-0,020225 

IhOoO 

0 11 

1.50 

3,398998 

0,145230 

-0.135782 

-0,129491 




m 

x= a. 

0050 



H 


u3 

K 

K 

a 

q tX,0) 




or 

oi 

fc 



10000.0 

0,75 

1.530119 

0,015450 

-0,011072 

lUUUO.O 

0,85 

1. 745333 

0.013178 

-0.008487 

louuu.o 

0,,9 5 

1,958616 

0,016494 

-0.011400 

11500,0 

O', 65 

1,342742 

0.010428 

-0,0^6710 

11500,0 

U . 7 0 

1,453813 

0,007308 

-0 ,003467 

1150'', 0 

0,80 

1.6/4846 

0.00*979 

-0.000831 

11500, 0 

0,«o 

1 ,893969 

y,u068/7 

-0,004332 

11500,0 

1,00 

2.110832 

0.U2U031 

-0.0149/B 

13000.0 

0.55 

1,142201 

0,009732 

-0.006544 

13000,0 

U.bO 

1,255935 

0,005120 

-.0.001833 

1 3000,0 

0 , / 0 

1.482893 

-0 , 0006b4 

0,004107 

13000.0 

0,7 5 

1, 595869 

-0,001222 

0.004924 

13000,0 

U.8 0 

1,708388 

-0,000191 

0,004072 

1 3000,0 

0,90 

1,931789 

0,007550 

-0,003233 

1 3000,0 

0.95 

2.042563 

0.014537 

-0,009958 

1 1 00 0.0 

0,50 

1,040903 

0,009550 

-0,006631 

14000,0 

0 • oO 

1,2/1181 

0,000034 

0.003082 

1 1000,0 

' 0 • / 0 

1.500960 

-0,004405 

0.007793 

3 4000.0 

0,8 0 

1,729182 

-0.002072 

0,005817 

14000,0 

0,90 

1,955160 

0.008395 

-0,004189 

1 oOUO.O 

0,40 

0,845718 

0,004794 

-0,002557 

3 8(!a0.0 

0,60 

1.325340 

-0,013001 

0,015658 

1 8‘'»U0 , 0 

0,80 

1,802410 

-0,000282 

0,003696 

1 yoyu , 0 

0.90 

2.036809 

0,022242 

-0.016245 


-0,002iB7 
"Ot unb228 
"tJ , UO 1^49 
0,001629 
0 . 00433 ^ 
u . uoa 8 y 4 
-0,00970? 
-0,002196 
0.002641 
0,OOB53B 
0,009386 
0, 0085ol 
0,001305 
-0,005399 
-0.002653 
0 . 00/102 
0,011866 
0,009941 
- 0,000022 
0.000418 
0, 018682 
0,000790 
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